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CHAPTER I 


INTRODUCTION ' 

Before the invention of the digital computer, elaborate and com- 
plicated numerical techniques for solving problems in science, mathe- 
matics, and engineering were only given secondary consideration. As 
the refinement of the digital computer progressed, its comprehensive 
usefulness became more apparent. Today, the employment of the digital 
computer is found in almost every discipline of science and engineer- 
ing. 

Mathematical Programming 

One area in which the digital computer has been of tremendous aid 
is in the solution of mathematical programming problems. The general 
mathematical programming problem may be stated as: 

T 

determine the n component vector x = (x , ■ x 2 , , . „ , x ) so that 

T 

the maximum (or minimum) of f (x ) 

T 

subject to g A (x ){<_, =, >_} c ± 

i - 1, 2, .... m (1-1) 

is obtained. Each relation in (1-1) is assumed to be algebraic in 

T 

nature. The relation, f(x ), is called the cost function, whose 
extremal with respect to the m constraints of the second relation is 
desired. If all the functions in (1-1) are linear and if the variables 
are not required to be Integral valued, then the above optimization 



problem is said to be a continuous linear programming problem (LP), 

The solution of the continuous linear programming problem may be 
accomplished with the aid of the simplex algorithm first introduced 
by George Dantzig. 1 Today the solution of continuous linear pro- 
gramming problems is treated extensively in many text books . 2 ’ 3 ' H ’ 5 

On the other hand, if any of the algebraic functions in (1-1) are 
nonlinear, then the problem is called a nonlinear programming problem 
(NLP). Up to now there has been no one algorithm developed that will 
solve all nonlinear programming problems. Generally the existence 
and uniqueness of a solution cannot even be assured without the cost 
function and the constraints possessing certain convexity and con- 
cavity properties. By placing various restrictions on the functions 
in (1-1) , there have been several algorithms developed for obtaining 
solutions. 6 In general, NLP algorithms are classed as either simplex 
in nature or as gradient in nature. 

Simplex Algorithms - Probably the first NLP algorithms developed 
were the separable programming algorithms. Problems for which they 
are applicable are of the following form: 


l g A * =, >}c , 


j=l 


'ij j ~ ~ i 


0 j = 1, ■*. . . , n 


maximize (or minimize) z - [ f.(x,) 

j=l 3 3 


(1-2) 


Here dynamic programming is not considered as a NLP algorithm 
but is considered as another branch of mathematical programming. 
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In order to apply separable programming both the constraints and the 
cost functions must be separable into functions of single variables. 

The mono-variable functions are then approximated over some finite 
interval by sequences of straight lines. Then a simplex algorithm is 
used to solve the approximate problem. The separable programming 
algorithms differ in the way the approximations are made and in the 
type of simplex algorithm necessary to solve the problem. 6 

Another simplex NLP algorithm is the quadratic programming of 
Wolfe. 11 It was especially developed to solve problems of the form: 

Ax = b 
x >_ 0 

T 

maximize (minimize) z = cx + x Dx (1-3) 

where A is an mx n matrix, c is an nx 1 matrix, and D is an n x n 
negative semidefinite matrix. In this case the constraints are linear 
and the cost function is quadratic and concave. The development of 
the algorithm for solving (1-3) depends heavily upon the Kuhn-Tucker 
conditions. 2 » 7 » 1 J 

Still another simplex type algorithm is the Hocking-Hartley con- 
vex programming technique^ 2 It is used to solve general NLP programming 

* 

problems with certain convexity and concavity conditions. This method 
is derived by approximating the cost function and the constraints by 
an infinite number of supporting hyperplanes „ Of course this produces 
a LP with an infinite number of rows. Then by using the duality prin- 
ciple of LP the problem is transformed into an infinite column problem 
which is amenable to solution by the simplex method. The convergence 
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properties of this algorithm are very reminiscent of the Newton-Raphson 
method for finding the roots of a polynomial, i. e. , whenever the 
algorithm converges, it usually converges very rapidly. 12 

Of course, there are many other simplex type NLP algorithms; in 
fact, there are several versions of those given above. However, for 
brevity only the more publicized algorithms and the basic thoughts 
behind them have been mentioned here. 

Gradient Algorithms — In contrast to the simples type algorithms 

i 

there exist the gradient algorithms. The premier algorithms of this 
type are the gradient projection method 13 » 14 , the generalized reduced 
gradient method (GRG) 18 , and the sequentially unconstrained minimi- 
zation technique (SUMT) 1 5 » 1 6 » 17 

The basic idea of the gradient projection method is to start with 

a feasible solution and move in the direction of the gradient of the 

cost function (for maximization problems) until the solution is found 

* 

or until the violation of a constraint is attempted. If the viola- 

I 

tion of a constraint is attempted, a direction is determined so that 
an increase in the cost function results and no violation of the 
constraints occurs. If no direction can be determined then the 
solution has been found. 

In the case of linear constraints this simply requires projecting 
the gradient Into the space defined by the Intersection of all con- 
straints which are equalities at the point under consideration. This 
is done by determining 

r = Pd (1-4) 

A feasible solution is any point where no constraint is violated. 
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where r is the directional vector which points in the direction to 
move, d is the gradient of the cost function, and P is a projection 
matrix. The projection matrix is determined as 

P = I - Q(Q T Q) -1 Q T (1-5) 

where Q is a matrix whose columns are the gradients ef the constraints 
which are strict equalities at the point of question. Of course if 
Q becomes square then P = 0. This does not indicate a solution but 
simply indicates that the feasible solution is located at a cojrner of 
the solution space. (For determining the projected gradient for this 
case, see Hadley [6], p. 167). 

Another so-called gradient NLP method is the GRG method mentioned 
previously. This technique is a natural extension of the reduced 
gradient method of Wolfe to include nonlinear constraints. The 
reduced gradient method was developed to determine relative extremals 
of 


maximize 

f(x) 




subject to 

Ax < 

b 


(1-6) 


X i^ 

0 

x 1,2,... ,n . 



It is assumed that any n-row submatrix of A has rank n. Next, A is 
partitioned into an nx n submatrix C and a submatrix D, and b is 
similarly partitioned into c and d. Then slack variables y and z are 
added so that the constraints in (1-6) become 

Cx + y = c (1-7) 

Dx + z = d . (1-8) 
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All the constraints which are equalities are included in the C matrix. 
The variables of z are considered as dependent and those of y as 
independent. From (1-6), (1-7), and (1-8) it is easily seen that 



Ax 

= 

- C -1 Ay 

(1-9) 


Az 

= 

DC -1 Ay 

(1-10) 


V y f(x) 

= 

- Vf (x)C -1 

(1-11) 

where Ax 

and Az represent 

the 

changes in the x's and y's. 

V f(x) 

y 

is called 
« 

the reduced gradient 

and 

Vf(x) is the gradient of 

the cost 

function. 

From (1-9), (1-10), 

and 

(1-11) a set of rules has 

been 

devised for determining the correct changes in the x's and y’ 

s so that 

an increase in f(x) is registered 

(For additional information 

see 


[30]). 

Somewhat different from the gradient projection and GRG methods 
is the SUMT. The problems amenable to this technique are those which 
can be cast into the following form: 

T 

minimize f (x ) 

T 

subject to g^(x ) _> 0 , i a 1, 2, q 

hj(x T ) = 0 , j ■ 1, 2 p. . (1-12) 

In applying SUMT the above constrained minimization problem is trans- 
formed and solved as a sequence of unconstrained minimization problems 
which in the limit converges to a solution. This is done by forming 
from the above cost function and constraints a penalty function of 
the following form: 
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P(x T ,R) = f(x T )+R l ^+4 ! ^2 (x T ) (1-13) 

1=1 g.(x T ) R j=l J 

where R is a weighting constant greater than 0. For some initial 
value of R the unconstrained penalty function, (1-13), is minimized 
by some unconstrained minimization technique. Then R is decreased by 
dividing it by some number greater than 1 and the process is repeated. 
As R -*■ 0 the unconstrained solution approaches a constrained solution. 
The physical effects of the two latter terms in (1-13) is to penalize 
a trial solution for getting too close to the boundary of the feasible 
region. 

There has been no attempt here to be all inclusive with respect 
to gradient algorithms. There are several other gradient algorithms 
that have been developed. However, the ones mentioned above are con- 
sidered by many as the most prominent and useful methods today. 

J ; r ‘ * 

Mathematical Programming in the Design of Control Systems 

Over the past ten years there has been a great thrust to use 
mathematical programming in the design of control systems. The major 
effort has been in the solution of optimal control problems, and the 
results in this area have been very fruitful — not only in the appli- 
cation. of mathematical programming but also in theoretical develop- 
ments.- In fact, it has been shown that the Kuhn-Tucker necessary 
conditions of mathematical programming and the maximum principle of 
optimal control can be derived from the same set of general optimiza- 
tion theorems 19 ,2 ° > 21 ’ 22 . As can be seen from the lengthy reference 
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list by Tabak 23 , much of the work has been directed toward the appli- 
cations of linear and quadratic programming.. Recently, uses of the 
SUMT and the GRG algorithms in the solution of optimal control 
problems have been made, 21+ » 25 »- 26 

On the other hand, the use of mathematical programming, in the 
classical design of control systems has been meager-^-particularly in 
the design of. compensators. from a frequency domain point of view. 

This is very -unfortunate because. most practical system designs even 
today are still by classical frequency domain approaches. . Further- 
more, these approaches are more artful than analytical. The. few 
techniques which -have been deyeloped can be classified as modern con- 
trol oriented- or strictly classical control oriented. This classifi- 
cation results, from the choices of the performance indices. Those 
methods in which system specifications are submerged in a cost 
functional are labeled as modern control approaches, while those 
methods which represent the system performance by classical standards 
such as gain margins, phase t margins , bandwidth, etc. , are termed as 
classical approaches. 

One of the first successful computerized compensator algorithms 
was developed by Coffey, 27 In his paper consideration is given to 
a system similar to that shown in Figure 1. In this figure j parame- 
ters of the system are sensed; each parameter is operated on by some 
compensation device; the results of these are summed and fed back. 
Figure 1 is considered typical of large aircraft or space vehicles. 
Each compensator is assumed in the following form : 
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(1-14) 


M 

e 

G e (s) = [ .\ a ei s 

i=l 


N 


i-1 


]/U+ I b. 


i-li 


i=2 


ei 


'j- 


where s. is a complex variable and M - 1 and N - 1 are the nume- 

e e vt’v 

til 

ra tor. and the denominator orders, respectively, of the e L compensator. 
The goal is to select the compensator coefficients so that the com- 
pensated open loop frequency response is a weighted least squares fit 
to a. desired open looj> frequency response (Of course, the open loop 
frequency response is obtained by calculating C(ju))/R(j<jj) when the 
feedback loop is broken at a). 


The weighted least-squares fit is obtained by minimizing the 
following cost function: 

J = | (y* - y*) T WW T (y - y) | (1-15) 

- A 

where y is a vector of the desired frequency response points, y is a 
vector of frequency response points, W is a diagonal weighting matrix, 
the asterisk (*) denotes conjugate, and the super T denotes transpose. 
For minimizing the cost function, J, with respect to the compensator 
coefficients a gradient search algorithm is chosen; and, then, the 
necessary gradient vector is calculated. 

Next, geometrical properties of the cost function are considered. 
It is demonstrated that even for relatively simple systems the cost 
function is geometrically complicated. From this it is seen that the 
cost functions can have relative extremals and unbounded solutions. 
Furthermore, the design of unstable compensators is possible. Even 
with the possibility of these difficulties, it is demonstrated that 
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this procedure can be utilized to design practical compensators. This 
is done by applying the technique tc a sixth order ballistic missile 
example. For this system two compensators are designed — .a pure gain 
and a fourth order over a sixth order. The pure gain compensator 
approximated th'e desired frequency response for low frequencies but 
was completely unsatisfactory for higher frequencies. In fact, for 
this compensator the closed loop system is unstable. On the other 
hand the higher order compensator exhibited very good properties -when 
compared to the desired frequency response. 

Coffey indicates that in some instances a judicious choide of 
the elements of the weighting matrix, W, is required before an 
acceptable design can be achieved. Thus, a computer program of this 
algorithm might require several runs — while juggling these elements 
between runs — before the proper values are conceived. Even with this 
disadvantage the algorithm is definitely superior to classical means. 

Another technique for computerized design of compensators for 
control systems has been presented by Page and Stear. 28,29 The thesis 
of this procedure is to vary the compensator coefficients until 
certain chosen frequency response specifications are satisfied. The 
procedure for attempting to do this is 

N . 

minimize F = £ K. (1 - S^/S. ) 2 (1-16) 

i=l 

& 

where N is the number of specifications considered, S_^ is the speci- 
fication as a function of the compensator coefficients, is the 

desired specification, and is a weighting constant. The constant 

■ sl d 

is chosen as positive, in general one, for <_ and as zero 
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a d 

for. S-£ > . This results in a satisfied specification being 

neglected. The goal is to drive F to zero. The reason for the choice 
of the above criterion function (1-16) is to try to. place the most 
emphasis on the specifications which have the greatest violations. 

In order to illustrate the given procedure Stear and Page pre- 
sent an example of the design of an autopilot for an aircraft. In 
accomplishing this design four unconstrained optimization procedures 
are used. Three are local search procedures, and one is a global 
search technique. As in the case of Coffey's cost function it is 
discovered that even for simple compensators the specification 
function (1-16) has relative extremals. From this it is deduced that 
the global search procedure is more applicable than the local search 
techniques if the starting compensator is strictly arbitrary. However, 
if a priori knowledge is used in picking the initial compensator this 
deduction is not necessarily true. 

Pitfalls of Previous Works on Computerized Compensator Design 
Procedures 

The two previously mentioned works on computerized compensator 
design procedures suffer from several drawbacks. First the procedure 
presented by Coffey is basically a frequency response shaping technique. 
In the design of compensators for most control systems, this is too 
rigorous; i. e., this requires the compensator to satisfy more con- 
straints than are necessary. Thus, the probability of all system 
specifications being satisfied is less. Another interesting fact is 
that in many instances the frequency responses of control systems 
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are not required, to match a desired frequency response--f requeqcy c to 
frequency — but'"are desired to have some general shape which . can be 
translated with respect to frequency. Even more conceivable is the 
desirability to -have several bands of the frequency response to be 
various distances from the -1 + jO point of the GH(jw) -plane and to 
have other bands of the frequency response constrained to be greater 
than or less than limitations with -respect to the origin of the 
GH(jw) -plane. Constraints such as these are not as strenuous as those 
requiring the frequency response to fit closely to some desired fre- 
quency response. 

A pitfall which is common to both the ? Coffey method and the 
Stear and Page method is the necessity of choosing some constants — in 
particular, the elements of the diagonal matrix, W, and the K^'s. It 
is obvious that in many situations a judicious choice of these must 
be made before any useful results will emerge. It was suggested by 
these authors that computer programs containing the algorithms may 
require several runs with various values of these constants before an 
acceptable design is achieved. However, this involves trial and 
error which was one of the justifications for going to a computerized 
procedure. 

Another drawback of the two algorithms presented is that some 
specifications may become worse while others become better. This 
immediately poses some serious questions, such as, what is a reason- 
able trade-off and where does it exist? If minimum standards of 
system performance have been set, it is very probable that nothing 
short of these are acceptable. In this case there is no trade-off. 
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On the other hand, it may be viewed that in practical designs it is 
not unusual to accept performances a little less than that desired. 

In instances such as this, performance tolerances must be set. 

Another shortcoming of the two methods is their failure to 
include inherent devices for maintaining compensator stability. If 
the designed compensator is unstable, then the stability criterion of 
the system changes completely. The result might be system instability 
which removes the compensator from the realm of a practical design. 

'• -V.. ’ . 

What is needed is an algorithm which tends to improve system specifi- 
cations at every iteration. Of course this might require the allowance 
of only incremental changes in the compensator coefficients. 

Another pitfall of the two previously mentioned works is the lack 
of any theoretical inclusions on compensator limitations. That is, 
none of the authors presented any theoretical developments showing 
what could be expected from their algorithms for a certain compensator 
order in a particular system. Thus, initially there is no way to know 
what minimum amount of compensation is necessary. In addition, these 
works presented no theory which indicates that the algorithms will 

produce a final compensator that is any better than the initial com- 

* 1 ; 

pensator . 

In essence, the techniques of Coffey and Stear and Page are 
"firsts" in the use of the computer for compensator designs, but they 
are somewhat limited. They do not present universal solutions in 
regard to computerized-compensation. It is the purpose of this dis- 
sertation to present the theory and a method of computer-aided 
compensator design that does not have the drawbacks of the previously 
presented techniques and is thus more universal. 
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CHAPTER II 


-FREQUENCY RESPONSE CONSIDERATIONS IN THE 
: ? DESIGN OF A CONTROL SYSTEM 

i 

Before the design of a system can be accomplished, the limitations 
or constraints and the desired performance of the system must be 
established. The measurement of the performance of the system is 
determined by comparing that obtained to that desired. Because of the 
limitations, in many instances, the desired performance cannot be 
achieved. In designing compensators for practical control systems 
there are, in general, two types of performance indices — time domain 
indices and frequency domain indices. Although it is quite obvious 
that these are related, no analytical means, up to now, have been 
devised for defining this relation except for the simplest control 
systems — less than third order. In practical designs the main limi- 
tations are system stability, nonlinearity, time variance, and 
sensitivity. Today many systems are designed by using linearized 
frozen time models and applying frequency domain concepts. 

Concept of Relative Stability 

In most practical systems stability is a major constraint. In 
fact, in most system designs a specified degree of stability is 
required. A specific degree of relative stability is required because 
of inaccuracies in the model of the system or in order to deter insta- 
bility if future parameter variations in the system plant result. 


Sometimes a certain amount of relative stability is desired to keep 
the system from resonating unnecessarily. 

In the past the degree of relative stability has been denoted 
by the classical gain (GM) and phase margins (PM) . However, in some 
instances these can be very misleading. For example, consider the 
hypothetical s-plane frequency response shown in Figure 2 which 
possesses acceptable classical stability margins (GM _> 2.0, PM >_ 30°) 
but which comes within some small distance of the -1 + jO point. Such 
a condition could represent a system which was very close to insta- 
bility. A better measurement of relative stability is defined as 
follows : 

A stability margin is defined as the magnitude of the 1 + GH(jw) 
frequency response at one of its minima relative to the origin 
of the 1 + GH(joj) plane. 

It is deemed by this author that by measuring stability in this 
fashion, a measure of the true relative stability of a system is 
achieved. Next, a system is said to be relative stable if the fre- 
quency response does not cross a designated closed contour located 
around the -1 + jO point. This closed contour around the -1 + jO 
point is called the margin of stability limit. 7 The shape and the 
size of this contour depend upon system specifications. Furthermore, 
there is nothing wrong with making the size and shape of the contour 
frequency dependent. (In doing this the designer would be indicating 
that the frequency response is to be shaped to some extent.) 

Relative Attenuation Concept 

Although relative stability plans a major role in compensator 
determination, there are several other factors which are considered. 
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Figure 2. A Hypothetical GH(jw) Frequency Response 


One of these is the attenuation of certain frequency 0 bands. The reason 
for frequency band attenuation is to discourage the control system 
from resonating at some natural frequency of the sys'tem. Of course if 
the system is linear and time-invariant this is not necessary. Un- 
fortunately, many practical systems do not fit into the linear, time- 
invariant category. 

Frequency band attenuation may be treated by requiring that all 
frequency points that are to be attenuated fall within a chosen con- 
tour around the origin in the GH(joo) plane. This contour is called 
the margin of attenuation limit . It then follows that: 

An attenuation margin is the magnitude of the GH(jm) 
frequency response at one of its maxima with respect to 
the origin of the GH(joj) plane . 7 

Other Frequency Response Concepts 

Relative stability and attenuation are considered as the most 
important frequency response design criteria. However, they do not 
yield acceptable designs in all instances. Sometimes it is necessary 
to employ proper phasing of certain frequencies. This is usually 
employed when it becomes difficult to determine a compensator ‘to 
attenuate certain natural frequencies of the system and in addition 

to satisfy other system requirements. The general idea' is to 

! : . . 

determine a compensator so that these frequencies are phased toward 

the right half of the GH(jw) plane.' This results in these frequencies 
being attenuated in the closed loop system; 1 

In some cases it is even necessary to place special- emphasis on 
certain points of the frequency response. In most instances these 
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points are closely related to dynamical responses of the control 
system. Examples of dynamical responses considered for a space craft 
are wind response and "engine^oub" response. In order that these 
responses possess acceptable characteristics it is usually rtecessary 
to require certain frequency response points to be placed in certain 
regions of the GH(jw) plane. 

Still another frequency response design concept is bandwidth. 
However, this can be handled by either the stability margin or the 
attenuation margins. For example, the maximum open loop bandwidth 
can be achieved by requiring a certain frequency and all frequencies 
above it to have a certain margin of attenuation limit. Similarly, 
closed loop bandwidth could be controlled by a combination of these. 

Problem Formulation 

Assuming that the desired frequency response characteristics have 
been determined so that if they are achieved the performance of the 
system will be acceptable, it must be decided how to determine a 
compensator for achieving these. The classical means of doing this 
is .by trial and error: however, a more efficient method would be an 

■ . ' j , : ■ ! 

iterative method that makes improvements upon the system's frequency 
response from iteration to iteration or indicates that no further 
improvement could be made. In fact, if a total of n critical frequency 

points have been chosen, then the problem may be formulated as the 

‘ « ", * ' 

following nonlinear programming problem: 

T 

Determine a vector x such that , , 
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( 2 - 1 ) 


gi (x T , o. i ) < b. 

*^ 1 * =» d ± 

, i — 1) n 

T 

In (2-1) x is a vector of the compensator coefficients; is a 
th 

function of the i frequency, to^, and the compensator coefficients. 

The functions , g^, i = 1, n, are chosen so as to represent the 

frequency response limitations and constraints which have been imposed. 
For example, gj, could be representative of a stability margin or an 
attenuation margin. The second relation in (2-1) takes into account 
any constraints that might be placed on the compensator coefficients. 

It may be necessary to constrain some of the coefficients if it is 
desired to keep; the d. c.. gain, G(jO), of the system constant or above 
or below a certain level. Also, it may be necessary to constrain 
certain compensator coefficients to insure the stability of the 
compensator or to take into account realizability conditions. 

The above formulated nonlinear programming problem differs from 
the classical nonlinear programming problem in the respect that it is 
strictly a constraint problem. 6 There is no cost function to maximize 
or minimize. However, this does not simplify matters. In fact, the 
above problem can be thought of as a normal nonlinear programming 
problem in which it is desired to find a solution which obtains a 
certain objective function value. In this case the objective function 
just, becomes a constraint., If the objective function is added to the 
constraint, list, then the result is a strict constraint problem as 
given above. The desired solution to this problem is a feasible 
solution which may or may not exist. 
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CHAPTER III 

COMPENSATOR LIMITATIONS 

At any iteration in solving the problem mentioned in Chapter II, 
there will result conditions of the form of (2-1) to be improved. 

(The number n can change from one iteration to another since the 
frequency response changes with respect to -the compensator.) The 
general idea is to change the compensator coefficients so that each 
constraint comes closer to being satisfied; The question then is, 
how many compensator coefficients are required to insure that some 
improvement on each constraint at a certain iteration can be made? 
This question is answered by the following definitions and theorems. 

Definition _1 

An optimal direction in the GH(jw) plane is any chosen 
direction in which it is desired to perturb a point on 
the frequency response. 

, J 

Optimal directions are illustrated in Figure 3 at points A, B, and 
C. The number of compensator coefficients sufficient to perturb n 
polar frequency response points in n optimal directions is given by 
the following theorem: 

t 

Theorem 1^ 

A sufficient condition to perturb n points on a polar frequency 
response curve in n optimal directions with a realizable compensator 
is that there be at least 2n independent compensator coefficients 
which are available to be varied. 
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Figure 3. A GH(jco) Frequency Response for Illustrating 
Optimal and Sub-optimal Directions 
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Proof : Let the open loop frequency response be denoted by 

T T 

G (jui, x ) where x is an m dimensional vector of the functionally 

independent compensator coefficients. Also, let the optimal 

* 

direction at a frequency a> k be denoted by d^ . Suppose there are n 
points on the frequency response which are to be moved in the n 
chosen directions, respectively. The change of the open loop transfer 
function at the k 1 "^ 1 frequency with respect to the i^ compensator 
coefficient is of the form 


3 < v«v * > 


3x . 


= C ki + j \i 


(3-1) 


where c^ and e^ are real constants. There are, for a particular 
frequency, m such partials as (3-1) and, if they were included as the 
components of a single vector, the result would be the complex 
gradient. It is well known that this points in the direction of the 
most rapid change. However, this is not the desired direction of 
movement. Essentially what is needed is a directional vector [w] in 
complex m-space whose dot product with the m dimensional gradient 

£ 

vector [c^ + je^.] will yield the desired directional derivative d^ , 
or in equation form (See [32]) 

d k = [c k + je k ] T [w] . (3-2) 

It should be obvious that the components of [w] are proportional to 
the amount that each compensator coefficient must be varied in order 

■k 

that movement in the d k direction can be accomplished. Thus if the 
compensator is to be realizable, [w] must be a real vector. 


23 


Letting 


d k " a k + ^ b k * 


(3-3) 


then (3-2) can be written by the following two real equations: 


and 


* T> 

\ = M 


\ ' M ■ 


(3-4a) 


(3-4b) 


Hence, for n points on the frequency response to be moved in n 
optimal directions there result 2n equations or 


* T 

= [Cj] [w] 

& T 

b 2 = t c 2^ l w ^ 


* 'p 

a = [c ] l [w] 
n 1 n J 1 


* T 

bj = [e-j] 1 [w] 


b 2 * = [e 2 ] T [w] 




(3-5) 


In matrix notation (3-5) becomes 


[w] 


(3-6) 


where the dimensions of 


* 

a 

• • • « 
L b* 


L e 


, and [w] are 


respectively 2n x 1, 2n x m, and m x 1. If 2n > m there will result 
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more equations than unknowns and possibly an incompatibility. 31 ' Hence 
there may not exist a vector [w] such that all equations can be satis- 
fied. This says there are not enough compensator coefficients avail- 
able.' On the r otlier hand if 2n £ m, there either results less equations 
than unknowns or the same equations as unknowns. For the first case 
there will exist an infinite number of vector [w] ' s and an infinite 
number of solutions to the equations. This indicates an excessive 
number of compensator coefficients. In the second case there will be 
a unique [w] and, thereby, a unique solution for the equations. This 
means that the exact number of compensator coefficients necessary is 
being employed.® 

The preceding proof has shown the sufficiency condition for mov- 
ing the frequency response in n optimal directions. Suppose, however, 
that it is desirable to use a compensator with a fewer number of 
coefficients than those needed to move in the optimal directions. 
Consider the following definition: 

Definition J2 

A sub-optimal direction is any direction within tt / 2 

radians of an optimal direction. 

An optimal direction is just a two-space vector; then, a sub-optimal 
direction is any two-space vector which has a positive dot product 
with an optimal direction. Thus, a sub-optimal direction is any 
vector which falls within a certain open half space, e.g., a sub- 
optimal direction to B in Figure 3 is any vector which points to the 
left of the line passing through B. 

If the optimal and sub-optimal directions for are respectively 
represented in 2-space by the following vectors : 
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(3-7) 



9 


and _ 

: O - • A " < a k ’ V * (3 " 8) 

then the sub-optimal direction would be any direction such that the 
dot product 

■ : • • ■+■••• ■* * 

d k • d k >0 (3-9) 

or 

Vk* + Vk* > 0 • (3 ‘ 10) 


Then the question is, how many compensator coefficients are necessary 
in order to assure that movement' in some sub-optimal direction can be 
achieved? The answer to this is stated and proved in the supervening 
theorem. v 


Theorem 2. 

In order to be assured of perturbing n points of an open loop 
frequency response in n sub-optimal directions, by varying the compen- 
sator coefficients, it is necessary that there be n independent 
compensator coefficients available for variance. 

til 

Proof : The components of the k C sub-optimal vector direction in 

til 

terms of the real and imaginary parts of the partials at the k 
frequency are given by 


m 


l c. , w 


i=l 


ki i 


(3-11) 
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(3-12) 


m 


\ e ki”i 


where c^ and e^, respectively, are the real and imaginary parts 
(evaluated at ui^) of the partial of the open loop transfer function 
with respect to the i tb compensator coefficient, and w^ is the i tb 
unknown constant which is to be determined so that (a , b, ) points in 

K K 

a sub-optimal direction. Substituting (3-11) and (3-12) into (3-10) 
results in 


m 

I 

i=l 


C ki w i a k 


m 

I 

i=l 


e ki W i b k 


> 0 


(3-13) 


or 


m 

l 

i=l 


l (c 


ki k 


+ e 


ki 


b k*> "i 


> 0 


(3-14) 


Remembering that there are n frequency points, n inequalities 
like (3-14) will result. Hence the following matrix inequality can 
be obtained: 


[c^a + e T b ] [w] > 0 . (3-15) 

The dimension of [c a + e b ] is n x m. In order to be assured that 
all n inequalities can be satisfied, it is necessary that there be at 
least the same number of unknowns as inequalities. Hence, this says 
there must be at least n independent compensator coefficients in order 
to be assured that n frequency points can be perturbed in the sub- 
optimal directions. 8 

The above two theorems place limitations on the overall compen- 
sator order. Thus for any algorithm to be assured of being able to 
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make the changes given in the theorems, the theorem must be sat- 
isfied. 


CHAPTER IV 


CONSTRAINT IMPROVEMENT ALGORITHM 

It is very desirable to have an algorithm which starts with 
some initial compensator and, then, in an iterative fashion produces 
an improved frequency response. This statement immediately suggests 
the question — what is an improved frequency response? This is 
answered by the following two definitions. 

Definition _3 

A total improved frequency response (TIFR) in an iterative 
scheme is one whose unsatisfied constraint values at a 
certain iteration are better than they were at the last 
iteration. 

Definition 4_ 

A sum improved frequency response (SIFR) in an iterative 
scheme is one whose sum of the differences in the unsatis- 
fied constraint values and their desired values is a positive 
value from one iteration to the next.* 

It is obvious that an algorithm which is capable of producing a TIFR 
is also capable of producing a SIFR; however, this statement is not 
reversible. A TIFR algorithm requires every constraint which is 
unsatisfied to be improved or bettered at every iteration, while a 
SIFR algorithm only necessitates a sum improvement, i. e. , the sum 


It is assumed, here, that all constraints in (2-1) have been 
represented in the £ form by multiplying >; constraints by -1 and 
changing = constraints to two inequality constraints (See Hadley 
[6]). No generality is lost by doing this. 
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increase must be better than the sum decrease. The goal is then to 
derive an algorithm which is compatible to both TIFR and SIFR. 

Thus, an algorithm is needed for solving a nonlinear program- 
ming problem of the following form: 

T 

Determine the vector x such that 

g^x 1 ) 2. b i 1 = 1* ••• , m . . (4-1) 

Again this is strictly a constraint problem. If .this problem has a 
solution, then it is a point in a solution space (Theoretically the 
solution space could be a single point). The functions in (4rl). are 
not assumed either concave or convex. What is desired is an iterative 
algorithm which, when started at some initial guess at the solution, 
will at each iteration produce an improved solution from the solution 
at the last iteration or will indicate that no further improvement 
can be made. An improved, solution is defined as one which brings the 
constraints closer to being satisfied. 

Constraint Improvement Algorithm Derivation 

T 

Suppose that some initial starting point, x , has been chosen. 

K 

Of the m constraints, let n be the number not satisfied by this point. 
The constraints not satisfied are defined as the active constraints , 
and those satisfied are called the inactive constraints . Let J, 
contain the index numbers of the active constraints, i. e. , 

J = {kj, k 2 , ...» k^}. Essentially what is desired is a directional 
vector, D, by which the vector x can'be changed, and it will be possi- 
ble to get an improved solution. This vector can be calculated as 

'< ■ ■ 

D = a l V 8k 1 + a 2 V 8 k 2 + ••• + • (4-2) 
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In (4-2) 


SI . . 


ro 


k^ * ^2 , • < • 5 k^ e J 


Vg k denotes the gradient of the constraint, corresponding to the 

T 

index evaluated at , and {a^} is a set of constants that are to be 
determined. An improved solution can be assured if the a's are 
determined so that i* 


D *' Vg k± >0 i = l,...,n . 


(4-3) 


In other words the maximum rate of increase of g^ at x^ is in 
direction of Vg ki> but an increase in g^can be registered by 
in the direction. of any vector which has a positive component 
direction of the gradient. In fact, suppose that a value for 
the dot products in (4-3) is chosen. Then (4-3) becomes 


the 

traveling 
in the 
each of 


D * V 8k x = c l 

D • Vg k = c 2 

• .- 

: D . * V 8k n = c n 


(4-4) 


where the vector c = (c , c 2> •••> c n ) contains the chosen dot 
product resultants. Substituting (4-2) into (4-4) results in the 
following set of linear equations, 


(V 8k 1 * V 8k^ a l + (V Skj ’ v Sk 2 > a 2 + ••• + V Sk n > a n = C 1 

( V 8k 2 * v gk^ a l + ( y§k 2 *-■ y 8k 2 > a 2 + •••,+ (Vg k2 * V §k n > a n = c 2 

(V gkn • Vg ki )a! + (Vg kn • Vgk 2 )a 2 + ... + (Vg kn * Vg kn )a n = c n . (4-5) 
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Using matrix notation (4-5) becomes 

[VG T VG]a = c , (4-6) 

T 

in which a = [a^ a^ ... a^ and VG is a matrix whose columns 
are composed of the gradients of the active constraints (The matrix 
[VG VG] is the Gramian matrix of the gradient vectors under con- 
sideration — see Hildebrand [31].). 

If the gradient vectors are linearly independent then 

T -1 

a = [VG VG] c . (4-7) 

Hence, this will yield a's for. a desired dot product between the 
directional vector D and each gradient of the active constraints. 9 
By moving in the direction of D then it is possible to improve the 
present solution.* 

Algorithm Summation 

Using the derivation and the preceding terminology, the con- 
straint improvement algorithm may be summarized as follows: 

*k+l = *£ + h[VG]a 

f 

T T til th 

in which x, and x are the solution points at the (k + 1) and 


In the above derivation the gradients were used. However, 
vectors in the directions of the gradients will suffice. In fact, it 
has been found in practice that unit vectors in the directions of the 
gradients are more suitable when the gradient magnitudes become 
disproportioned . The main advantage is a greater convergence 
rate. 
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iterations respectively, [VG] is a matrix whose columns are composed 

T 

of the gradients of the active constraints evaluated at x^ , 

T “1 

a = [VG VG] c 

• r > t 

where c is a column matrix of positive constants, and h is a positive 
constant. 

The choice of h (the step size constant) determines how much or 
whether any improvement in the constraints is made. In a compensator 
design program h also determines whether the program is a TIFR or SIFR 
algorithm. As a general rule small positive values of h produce a 
TIFR and larger values of h produce a SIFR. Of course there is a max- 
imum limit on h for producing a SIFR, i. e. , values of h above the 
maximum do not produce either a TIFR or a SIFR. On the other hand, 
negative values of h are out of the question since they tend to 
decrease the constraints — making them even worse. 

In addition to choosing h, a choice of the components of the c 
vector must be made. As has been pointed out previously, the com- 
ponents of c are the dot products of the directional vector, D, and 
the gradients of the active constraints. Thus by properly choosing the 
c's the amount of increase in some of the constraints can be, to some 
extent, controlled. In other words by judicious choice of the c's some 
constraints can be weighted more heavily than others. However, the 
actual amount of change in a constraint is related to h and the con- 
straints' partial derivatives. In practice it has been found that 
when using unit vectors in the directions of the gradients of the con- 
straints a good choice of the elements of the c vector is l's. This 
choice gives the best convergence rate. 
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■•On the other, hand, there is nothing wrong withjjmaking the c's 
dependent upon the constraint values, e. g., by letting a c decrease 
as its .constraint comes closer to being satisfied. t .However , as a c 
approaches zero the algorithm would tend to determine, a direction that 
was parallel to the boundary of the feasible region. Hence, the proba- 
bility of the constraint corresponding to this c becoming inactive 
decreases. Nevertheless, it has been discovered that in many instances 
that by holding the c's at respectable positive levels many of the 
constraints are driven to inactivity and they do not return to activity 
again. In this case the order of the matrix whose inverse is required 
can be reduced, whereas, if all constraints always linger in activity 
the order can increase if other constraints become active on higher 
iterations. 

Algorithm Limitations and Termination 

Next, attention is focused on algorithm termination. There are 
three conditions in which the algorithm will terminate. These are 

1. All constraints are inactive. 

2. One of the gradients of one of the constraints becomes zero. 

3. The gradients of the active constraints become linearly 
dependent. 

The first of these simply indicates that a solution has been obtained. 
The second and third represent relative extremal solutions. In fact, 
the second one shows that the solution point is a relative extremal of 
one of the constraints. On the contrary, the third termination con- 
dition indicates that at least one of the constraint gradients is a 
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linear c'ombinat'iori of the others' gradients or there are more active 
constraints ‘than- "there are variables (This could represent an incom- 
patibility condition.). Whenever 2 or 3 occurs either the solution 
obtained will have to be accepted or a new starting point will have 
to be chosen and the algorithm reinitiated. 


CHAPTER V 


GENERALIZED PARTIAL CALCULATIONS 

''•‘In essence, the goal of the designer is to pull and push various 
points on the frequency response until system specifications have been 
met or until no further’ improvements can be accomplished by the present 
compensator. In general, this can be accomplished by pushing and 
pulling the various points with respect to other points in the complex 
GH(jio) plane. For example, relative stability can be obtained by push- 
ing the points of the stability margins away from the -1 + jO point. 

On the other hand, the attenuation margins can be improved by pulling 
these points toward the origin. Similarly, proper phasing could be 
achieved by attempting to pull or push these points with respect to 
real axis points. Of course, in some specialized cases it may even be 
advantageous to pull or push a point with respect to more than one 
point. Regardless of whether a point is to be pushed or pulled it is 
necessary to know how these points change with respect to other points 
in the GH(joi) plane. This is especially true if the algorithm in 
Chapter IV is to be used in perturbing these points. 

A point can be pushed or pulled with respect to another point, 

- K, in the complex GH(joj) plane by varying the distance squared, 
d(m), between the point and - K. In order to determine how this dis- 
tance changes with respect to the compensator coefficients, con- 
sideration is given to the general feedback control system shown in 


Figure 1. The open loop frequency response of this system is deters 
mined by breaking the feedback loop at o and then calculating 


GH(joj) = C(jw)/R(jw) . . (5-1) 

Furthermore, to generalize even further in Figure 1 each channel's 

compensator is assumed to be made up of a product of sub-compensators, 
til 

i.e., the k channel's compensator is given as 

"k 

G,(s) = n G (s) , (5-2) 

i=l 

th * 

where n^ is the number of sub-compensators in the k channel. The 
uncompensated open loop state frequency response of the k^ channel 
with all channels opened is defined as 



a k (w) + jb k (w) 


(5-3) 


where a^ is the real part and b k is the imaginary part. 

From the above equations and statements it then follows that 


d(co) 


K + l {[a k (u>) + jb k (m)][ n G ki (jm)]} 
k=l i=l 


2 


(5-4) 


By assuming each sub-compensator to be a general rational function of 
the following form 


G 

qp 


(s) = 


n n-1 . 

n n-1 

m . m-1 . 

y s + y , s + . . . 
y m ■'m-1 


+ x 


+ Y 



(5-5) 


This is called the factored form of a compensator. 
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it becomes necessary to derive only how d(co) changes with respect to 
the coefficients of this general compensator, because the change in d(u>) 
with respect to any compensators’ coefficients will assume the same 
general form, only differing by the orders, n and m, and the numerical 
values of the x’s and y's. Since G^p(jw) is completely independent of 
all the other compensators, then it may be isolated from the others in 
(5-4). This is easily done by letting 


j’ \ 

A + jB = K + l {[a (u) + jb, (w)][ n G, (jw) ] } 
• k-i- i=l 1 • 

k^q 

and n 

q 

c + jd = [a (tu) + jb (w)] II G (jcu) 

q q q 1 

i^P 

Using (5-6) and (5-7), (5-4) is rewritten as 


d(aj) = 


A + JB + (c + 


j d)G qp (ja> ) 


(5-6) 


(5-7) 


(5-8) 


Substituting (5-5) into (5-8) and carrying out the necessary manip- 
ulations (5-8) evolves as 


( J Q C i X i + A E 2j y 2j " B J 0 E 2j+1 y 2j+l } + 


k 2 


d(oi) 


( J Q D i X i + A l E 2j+1 y 2j+l + B .| 0 E 2j y 2j 
% E 2j y 2j ) + ( J 0 E 2j+l y 2j+l ) 


(5-9) 
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where. -in- (5—9) .*a.y k = m/2 and ; p •■>m/2 - •!' • 
if Vi'-'.- m is even 

or' : =s.~- :• ; k = (m-i )/2 and p •=< (m-l )/2 

if ■ '• • \ -ji ; m is odd; ?'.*••« - • ; -• » yy>.- 

the C ' s i' D' s-,' and' E' s are defined by- the following sets': hit 

{Co> Cj, C 2 , C 3 , C 4 , C 5 , ...} = {c, -du, -c<o 2 , d<o 3 , co) 4 , -dui 5 , 

{Dq, Di, D 2 , D 3 , D 4 , D 5 , ...} = {d, co), -du> 2 , -cio 3 , dw 4 , cto 5 , . 

{Eq, Ei, E 2 , E 3 /E 4 , E 5 , ... } = {1, to, —to 2 , , — (o 3 , <o\ to 5 , ...}. 

(5-10a, 


Next, letting 


FNl - 

n 

1 

i=0 

k 

C i x i + A l 
j=o 

P 

E 2j y 2j " B E 2j+1 y 2j+l 

(5 


n 

P 

k 


FN2 = 

l 

i=0 

D x + A l 
1 1 • j=o 

E 2j+1 y 2j+l * B jl 0 E 2j y 2j 

(5 

FD1 = 

k 

l 

j=0 

E 2j y 2j ' 


(5 

FD2 = 

P 

I 

J-o 

E 2j+1 y 2j+l 


(5 

FD 

2 2 
(FD1) + (FD2) 


(5 

FN 

2 2 
(FNl) + .(FN2) 

1 1 

(5 


then 

9d(a) ) 2 [FD (A • FN1 + B • FN2) - FN • FDl]E q 

8 x (FD) 2 

q 


...} 

..} 

b,c) 

* 

-ID 

- 12 ) 

-13) 

-14) 

-15) 

-16) 

-17) 
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and 


where in (5-9) - ir > k = m/2 
if ' wi. m is even 


■p =,m /2 — - 1 * • 


or ' r: k = (m-l )/2 and p = '■ (m-l )/2 ? 

if ■ , .i 4 m is odd; . *- ••• •. • ! ? ,-v. : : 

the C's, D’s,. and' E's are defined by the following sets: ; 


{Cq , Cj, C 2 , C 3 , C 4 , C 5 , ...} = {c, -do), -cw 2 , do> 3 , co) 4 , -do) 5 , 

{Dq, Dj, D 2 , D 3 , D 4 , D 5 , ...} = (d, cu>, -du) 2 , -co) 3 , d(o 4 , co> 5 , . 

• ♦ . . . t 

{Eq, Ei, E 2 , E 3 , E 4 , E 5 , ...} = (1, u>, -oj f 2 ., -u 3 , w 4 , w 5 , ...}. 

(5-10a, 


Next, letting 


m = j V, + A t E y - E j E 


i =0 


n 


j =0 

P 


j =0 


J 2 j +1 y 2j+l 


FN2 = 


Jo °i X i + A .Jo E2J+1 72:5+1 + B Jo ^ 


2j 


(5 


(5 


m = Jo ^ ^ 


(5 


FD2 ~ J 0 E 2 j +1 y 2 j+l 

2 2 

FD = (FD1) + (FD2) 

2 2 

FN = (FN1) + (FN2) 


(5 

(5 

(5 


then 


gd^j 2[FD (A • FN1 + B • FN2) - FN • FDl]E q 


9x 


(FD) ‘ 


(5 


...} 

b,c) 

- 11 ) 

- 12 ) 

-13) 

-14) 

-15) 

-16) 

-17) 
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for q even or 


3d(m) 2[FD(-B • FNl + A • FN2) - FN. */FD2]E_ 

9x q (FD) 2 

for q odd and 


(5-18) 


3d(w) 2 [FNl ♦ c q + FN2 » D q ] 

3y q FD 


(5-19) 


for q even or odd. 

By programming "equations (5-4), (5-6), (5-10a,b,c), and (5-11) - 
(5-19) on the digital computer the partials of d(oi) with respect to 
the coefficients G (s) can be obtained. 9 ’ 10 

qp 

The above derivation provides the key for determining how any 
sub-compensator affects d(m) in a first order sense. With a complete 
comprehension of this derivation it becomes clearly apparent how to 
proceed either from channel to channel or from sub-compensator to sub- 
compensator in order to determine the necessary partial derivatives 
for a particular frequency point. Of course, this process must be 
completely repeated for each individual frequency point. Once the 
gradient vectors of each chosen frequency point are determined, then 
the calculation of the directional vector is accomplished as described 
in Chapter IV. 
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CHAPTER VI 


COMPENSATOR' IMPROVEMENT PROGRAM 

The preceding ideas were programmed in a digital computer program 
called CIP (Compensator Improvement Program). A complete fortran 
version of this program is contained in the Appendix. The general 
iterating procedure employed by CIP is as follows: 

1. Using the compensator at hand, the program calculates the 
critical points, i. e., stability margins, attenuation margins 
and other points of interest. 

2. If this is the first iteration a preselected step size is 
chosen. Otherwise, a step size is selected according to one 
of two criteria. 

3. Next the active constraints are separated from the inactive 
constraints. 

4. After this, unit vectors in the direction of the gradients 
with respect to the variable compensator coefficients are 
obtained (The numerator partials are listed first). 

5. Then using a chosen dot product vector the directional vector 
is determined (For the normalized gradient vectors calculated 
in 4, a suitable dot product vector has been found to be a 
vector whose components are l's). 

•k 

The other points of interest are frequency response points on 
which special attention is to be placed, for example, points to be 
properly phased, certain gain or phase margins, etc. 
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6. Finally, the directional vector is normalized with: respect 
to its magnitude; the compensator coefficients are changed 
according to the normalized directional vector and the step 
size; then, the complete process is repeated. • 

In order to initiate the program, an . input of discrete open loop 
frequency responses in the form of frequency and real and imaginary 
parts are required. Allowances are made for five channels of such 
information with a maximum of 999 points for each channel. This means 
that in Step 1 the actual critical points of the frequency response 
are not located — only approximate values are found. However, exper- 
ience has shown that the approximate values suffice. 

In order to determine better approximations to the critical 
points the input would require- open loop transfer functions (Equation 
5-3) for each channel. The more accurate approximations of the 
critical points could be found by finding the real roots of equations 
of degree 2n, where n is the total number of the open loop system (See 
5-1).* For systems above tenth order this is completely impractical 
due to the amount of computation time necessary to perform this task. 
Furthermore, in many practical situations an experimental discrete 
frequency response is the best information available for describing 
the system. In other words an experimental frequency response is 
obtained, and using this data a transfer function of the system is 
approximated. 

Also, some initial compensator for each allowable channel is 
required. The amount of initial compensation must be enough to 

*In this discussion it is assummed that due to round-off error 
a computer is not capable of getting exact solutions of non-integer 
problems. 
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stabilize the system.* If the system is open loop stable then each 
initial compensator can be chosen as an equivalent 1 compensator, i.e., 
the numerator and denominator factors are chosen to be the same. The 
compensators may be either in a factored or unfactored form (It is 
apparent that the unfactored form is just a special case of the 
factored form). 

In Step 2 the proper step size is chosen. In the CIP one of two 
procedures for selecting the step size is employed. These are ‘ 

a. Require the betterment of all active constraints from the 
last iteration. 

b. Require the sum of the differences of all active constraint 
values and their desired values to increase from the last 
iteration (For this sum all active constraints of the <_ 
form have been changed to the >_ form by multiplying by -1). 

Procedure a indicates the program is to be used in the TIFR phase, 
while procedure b designates the program as SIFR. The choice of the 
criteria used is left to the designer. If the one chosen is satisfied, 
the present step size is doubled, provided that the doubling process 
does not exceed some preselected maximum step size value.** Otherwise, 
the maximum step size value is utilized. Regardless of which of these 
occurs the program continues to the next iteration. On the other 
hand, if the continuance criterion is not satisfied then the step size 
is halved and the present iteration is repeated if the step size is 

* 

If the system is not stable then relative stability has no 
meaning — although relative instability might. 

•ff j|$f 

The main reason for limiting the step size is to keep the com- 
pensator from becoming unstable on a single iteration. 
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greater than some chosen minimum step size. When the step size 

* 

becomes less than the minimum value the program is terminated. 

Steps 3, 4, and 5 are simply operations necessary for employ- 
ment of the constraint improvement algorithm of Chapter V, whereas, 
in Step 6, the compensator coefficients are actually changed. In 
Step 5 the reason for reducing the directional vector to a unit 
vector is so that the step size actually designates the overall change 
in the compensator coefficients. Otherwise this would not be the 
case. 

The output of the CIP can be controlled to occur at every 
iteration or at set increments, i. e. , a set number of iterations 
can be skipped between outputs. At any iteration at which an output 
occurs the following information is printed by the CIP: 

1. Iteration number 

2. Constraint values 

3. Frequencies where constraints occur 

4. Desired constraint values 

5. Type of constraints 

6. Directional vector at the last iteration 

7. Compensators at the present iteration 

In 5 the type of constraints denotes whether it is a phase margin, a 
gain margin, a stability margin, or an attenuation margin, and the 
symbols used to denote these are respectively P, G, S, and A. 

In the program stability margins are the main vehicles for 
determining the relative stability of the system. The concepts of 

■k 

The program, also, has a maximum iteration termination condition. 
Since this has no effect on convergence, it was not included. 
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classical phase and gain margins have been included in the program 
because in some special cases these can be used to control proper 
phasing and various dynamical responses of the system. Furthermore, 
it should be pointed out that the measurement of these concepts is 
carried out exactly as stability margins, i. e, , distances from the 
-1 + jO point. Of course there is a one-to-one correspondence 
between this measuring method and the normal methods of measuring 
phase and gain margins. 


CHAPTER VII 


LARGE SYSTEM EXAMPLES 

In order to illustrate the practical usefulness of CIP, the 
improvements of the compensators for large systems are presented. 

This is done by way of two examples; the first example is a single 
channel system, while the second example is a dual channel system. 

The two systems are not the same, although they are very similar. 

Single Channel Example 

In this example the system under consideration is similar to 
that shown in Figure 1, but only one channel is fed back. The 
system's dynamics, 0i(s)/R(s), are described by the gain vs frequency 
and the. phase vs frequency plots shown in Figures 4 and 5. This sys- 
tem is a model of the Saturn V/Sl-C Dry Work Shop at a flight time of 
80 seconds. By an inspection of these frequency response plots it 
is revealed that this system has several poles near the j co-axis. 

This deduction is based on the spike shaped gain response and, the 
almost discontinuous changes in the phase response. These poles near 
the jco-axis are due to various sloshing and bending modes of the 
vehicle. 

This vehicle is inherently open loop unstable. Thus, it is 
necessary to use a control scheme, such as depicted by Figure 1, to 
stabilize it. Also, unity feedback with a pure gain compensator is 
not sufficient to stabilize the system. A compensator with unity 
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AMPLITUDE 



Frequency in CPS 


Figure 4. Gain vs Frequency for Uncompensated System 
of Single Channel Example 
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PHASE ANGLE IN DEGREES 



Frequency in CPS 


Figure 5. Phase vs Frequency for Uncompensated System 
of Single Channel Example 
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feedback which is capable of stabilizing the system is 


G (s) 
c 


0.9 


1.0 + 11.79440s +'28. 59200s 2 . 
1.0 + 21.56500s + 6.05650s 2 

1000.0. + 19.08700s 
1000.0 + 330.35200s 


100.0 4- 6.05720s + 7.56640s 2 
100.0 + 10.06500s + 6.32880s 


+ 3.73500s 2 
+ 19.02000s 2 


(El-1) 


2 


The GH(joj) compensated frequency response is shown in Figure 6. In- 
cluding the compensator, this frequency response represents a 29th 
over a 35th order system. 

In the design of the preceding compensator several physical 
limitations and constraints were considered — other than just stability 
of the system (In fact, stabilization of the system can be easily 
accomplished by a simple lead network with a reduced d. c. gain). 

Some of these are 

1. From past history it is known that compensators with very 
small d. c. gains produce poor wind responses. An acceptable 
value of d. c. gain is 0.9. 

2. On the GH(jw) frequency response the first negative real axis 
crossing with respect to increasing frequency is called the 
aerodynamical gain margin. Experimentation has shown that 
the major effect of an "engine-out" is a reduction of this 
margin. A safe crossing point is considered as -2 or less 
(or a frequency response magnitude greater than 2). 

3. For a small band of frequencies around 1.199 Hz the frequency 
response is dominated by the first bending mode. It is 
desirable to attenuate this band of frequencies. However, 

to even approach other system requirements and perform this 
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attenuation has been practically impossible. It has been 
found that the same effect results if this band of frequen- 
cies is phased in the right half of the GH(jw) plane. Due 
to the fact that the frequency of this mode is not known 
exactly, it is necessary to require larger phase margins for 
this mode than normally required. Acceptable margins are a 
lead phase margin of about 55° and a lag phase margin of 
about 90° (The reason for the difference is that in most 
physical systems phase lag is more probable to occur than 
phase lead). 

4. - For frequencies greater than 2.1 Hz the GH(jw) frequency 

response is dominated by the higher order bending modes. The 
control system can be deterred from resonating at any of 
these' higher modes by attenuating to a certain degree all 
frequencies above 2.1 Hz. These frequencies are considered 
satisfactorily attenuated if the magnitude of the GH(ju) 

\ 

frequency response is less than 0.25 for f > 2.1 Hz. 

i 

5. ' Besides the above frequency response requirements, it is 

desirable for all stability margins to be 0.5 or greater 
(Notice that in terms of classical stability margins this is 
approximately equivalent to having phase margins of 30° and 
gain margins of 2 or better) . 

By an observation of Figure 6 it becomes evident that all of the above 
specifications are not met. This becomes even more obvious after an 
inspection of Table 1. In this table the first margin is the aero- 
dynamical gain margin and the next two margins are the lead and lag 
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* 

phase margins of the 1st bending mode, respectively. The remaining 
margins listed under attenuated frequency information are stability 
margins as defined in this paper, and of course the attenuated infor- 
mation. is representative of the attenuation margins above f = 2.1 Hz. 

In the CIP program the following specifications were made: 

1. Determine the aerodynamical gain margin and improve it if 
it is less than 2. In order to improve any point it is 

i . • 

necessary to specify what point or points in the complex 
GH(jm) plane this point is to be pulled or pushed with 
j ^respect to. For this example it is chosen to push this 
point with respect to the -1 + jO point. 

2. Determine the lead and lag phase margins of the first bend- 
ing mode and improve either or both if they fall below 0.9 
and 1.3, respectively. To improve these it is chosen to 
push them from the -1 + jO point. 

3. Detect all stability margins and increase those less than 
0.505. Again the -1 + jO point is chosen as a pushing point. 

4. Detect all attenuation margins for f > 2.1 Hz. and decrease 

( 

all of' those greater than 0.25. For these margins the origin of 
the GH(jrn) plane is chosen as a pulling point. 


* 

The measurements of these stability margins are made in the 
same manner as stability margins defined in Chapter II, i.e. , the 
distance from the -1 + jO point. Measuring gain margins in this way 
is quite natural. However, measuring phase margins in this way is not 
as straight forward, even though there is a one to one correspondence. 
The equations relating the two are: d = 2 sin 0/2 and 6=2 arcsin d/2, 

where d is the distance from the -1 + jO point and 0 is the phase 
margin.- Of course d is limited to the closed interval [0,2]. 
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The continuance criterion chosen was b of Chapter VI. With these 
insertions and the necessary frequency response information in CIP, 
the following compensator was obtained after 2000 iterations (or 
approximately 30 minutes on a UNIVAC 1106): 

, s o 1.0 + 74.40524s + 107.13383s 2 100.0 + 7.29719s + 8.68710s 2 

V S; 1.0 + 124.68711s + 16.85849s 2 100.0 + 11.98668s + 9.15484s 2 

1000.0 + 12.10541s + 3.11162s 2 

. (El-2) 

1000.0 + 219.54201s + 20. 42297s‘ !: 

A tableau of the pertinent information at iteration 2000 is shown in 
Table 2. From this tableau it is seen that most margins are, for 
practical purposes, satisfied. The reason that several of the margins 
have values that are only approximately equal to the desired values is 
that, in most instances, after a margin becomes inactive it has a 
tendency to oscillate between activity and inactivity on higher itera- 
tions. However, by establishing an upper limit on the step size from 
iteration to iteration these constraints are coerced to remain in a 
vicinity of their desired values (For this example the maximum step 
size was chosen as 0.1 for the first 1000 iterations; then, to speed 
up convergence it was changed to 0.2 for the next 1000 iterations). 

The three smallest stability margins do not belong in the same 
category as those mentioned above because at no time were they inactive. 
Since program termination was maximum iterations, additional improve- 
ments in these constraints is quite conceivable. Nevertheless, the 
convergence curve shown in Figure 7 indicates many more iterations will 
be required before any appreciable change in the smallest of these 
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RELATIVE STABILITY INFORMATION 
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Figure 7. Convergence Curve for Single Channel Example 
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margins is recorded. With an occurrence such as this the designer is 
left with three alternatives: 

1. Accept the present design. 

2. Pay the toll of additional computer time and attempt 
additional iterations. 

3. Change some of the desired constraints and continue the 
program. 

From experience it has been found that small changes in the desired 

* 


margins can result 
cussion the GH(jw) 
practical purposes 


in marked effects. As for the case under dis- 
frequency response in Figure 8 reveals that for 
the compensator for iteration 2000 is satisfactory. 


10 


Dual Channel Example 


Again reference is made to Figure 1, except in this case it is 
assumed that j = 2, i. e., two channels are fed back. The uncompen- 
sated open loop system is described by the gain and phase frequency 
responses shown in Figures 9, 10, 11, and 12. Figures 9 and 10 
represent the gain and phase plots of 0!(s)/R(s), while Figures 11 
and 12 are the gain and phase plots of 02(s)/R(s). This system is 
typical of the Saturn V/Sl-C Sky Lab at a flight time of 105 seconds. 


* 

It should be noted that at the end of iteration 2000 the CIP was 
slightly modified so that a better calculation of the first negative 
real axis crossing frequency was obtained. After this, additional 
iterations were attempted and in less than 50 iterations the smallest 
stability margin was increased from 0.46513 to 0.48177. In another 
instance the compensator whose* smallest stability margin was 0.48177 
was used as the starting compensator in another run in which the rela- 
tive stability requirements were lowerecj to 0.49 while the other 
system requirements were the same as previously stated. In less than 
50 iterations all system requirements were completely satisfied. 
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90 ° 



270 ° 

Figure 8. GH(jaj) Compensated Frequency Response at 

Iteration 2000 for the Single Channel Example 
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Uncompensated System of Dual Channel Example 
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F^CHCNCr IN M7 


Figure 10. Phase vs Frequency for Channel 1 of Uncom- 
pensated System of Dual Channel Example 
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Figure 11. Gain vs Frequency for Channel 2 of Uncom- 
pensated System of Dual Channel Example 
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10. c 

10.0 + 2.85080s 


(E2-2) 


With these compensators Inserted In the system the compensated open 
loop GH(ja)) frequency response, C(jw)/R(jw), with the loop broken at 
a is that shown in Figure 13. 

It is desired to make several improvements in this frequency 
respons^. These conditional improvements are 

• j • • ' 

1. Keep the. aerodynamical gain margin at 4.37 or greater. 

2. Increase all stability margins of 0.49 or less. 

3. Maintain the lead and lag phase margins of the first 
bending mode at 55° and 90° or better. 

4. Decrease all attenuation margins occurring at frequencies 
above 2.0 hz when 0.2 or greater. 
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100 ° 270 ® 90 ® 280 ® 80 ° 


Figure 13. Initial GH(jw) Compensated Frequency Response 
for the Dual Channel Example 








In order to make these improvements the following specifications are 
made in CIP: 

1. Whenever the aerodynamical gain margin is 4.8 or less it is 
pulled with respect to the -7 -j3 point and pushed with 
respect to the -1 + jO point. 

2. All stability margin points less than 0.49 are pushed with 
respect to the -1 + jO. 

3. The lead and lag phase margins are pulled with respect to 
the 1 + jO point when less than 0.9 and 1.3 respectively. 
Also, the attenuation margins occurring at frequencies 
between these two are decreased by pulling with respect to 
the origin of the GH(jm) plane if they are greater than 9.0. 

4. The attenuation margins above 2.0 hz are decreased by pulling 
them with respect to the origin. 

With these specifications, 357 frequency response points for each 
channel, and the initial compensators, {E2-1) and (E2-2) , in the CIP, 
the following compensators were obtained after 200 iterations or about 
10 minutes on a Univac 1106 : 

„ , * . 1000.0 + 7.07293s + 7.02583s 2 100.0 

l? ( 3 ) ~ - ■ 1 ■ ■ — " ■ 1 — - - ■ ■ - ■ ■■■ 1 ^ ' 

2 100.0 + 1.21230s 100.0 + 10.48567s 

10.0 + 3.43938s 0.1 + 1.21370s 

10.0 + 1.14372s 0.1 + 2.51497s 

1.0 

10.0 + 1.14372s 


100.0 

100.0 + 9.34985s 


(E2-3) 
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1000.0 + 6.74527s +4. 53868s 2 100.0 + 4.57638s 

100.0 + 0.0s 100.0 + 1.26840s 


10.0 10.0 

100.0 + 6.96980s 10.0 + 1.44585s 


10.0 

10.0 + 1.44585s 


(E2-4) 


An evaluation of the amount of improvement can be made by comparing 
the initial tableau, Table 3, of important information to the final 
tableau. Table 4. As in the last example the first margin is the 
aerodynamical gain margin, and the next two margins are the lead and 
lag phase margins of the first bending mode respectively. The remain- 
ing margins under relative stability information are listed as stability 
margins. The margins under the attenuated frequency information are 
the attenuation margins above 1.2 hz. The desired margins' values are 
listed in the right hand column. 

Taking into account the desired improvements it is seen that 
significant improvement has been made. Furthermore, this is reinforced 
by comparing the initial compensated frequency response. Figure 13, to 
the compensated frequency response at iteration 200, Figure 14. The 
termination reason was maximum iterations; thus, as in the first 
example the designer is left with the same three alternatives. From 
the convergence curve shown in Figure 15, it appears that several 
additional iterations may have to be attempted before any significant 
improvement in the smallest stability margin is observed. The impor- 
tance of this example is the significant improvement over the initial 
frequency response. 
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RELATIVE STABILITY INFORMATION 
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Table 3. Initial Table of Values for Dual Channel Example 
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Table 4. Iteration 200 Tableau of Values for Dual Channel Example 






Figure 15. Convergence Curve of Dual Channel Example 
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Additional Analysis of Results and Comments 

The results obtained from the two examples clearly indicate that 
the CIP can be a valuable design aid. It must be pointed out that as 
the name, Compensator Improvement Program, imples the program is a 
design aid, not a design technique. That is, the program does not 
decide the order, the type, or the number of compensators necessary. 
All of this requires good engineering judgement before the running of 
the program is attempted. 

As the two examples exemplified the solution cannot be worse than 
the original compensator if the specifications on the input are made 
properly. In regard to stability margins and attenuation margins this 
simply requires pushing and pulling these, respectively, with respect 
to the -1 + jO and 0 + jO points. By doing this, these can always be 
bettered, except when they proceed from activity to inactivity. How- 
ever, the amount of slippage in going from inactivity to activity can 
be minimized by choosing a reasonable maximum step size such as 0.1 or 
less of the smallest compensator coefficient. As long as a margin 
stays in a vicinity of the desired value it is acceptable. 

The specifications for insuring the improvement in gain and 
phase margins are not always as simple as those for stability margins 
and attenuation margins. In fact, in many instances it is necessary 
to push and pull these with respect to two points in the complex 
plane. This is especially true if the acute angle between the tangent 
to the GH(jco) frequency response where these occur and either the tan- 
gent to the unit circle or real azis is very small. Both of these 
cases are illustrated in Figure 16 where tangents to some hypothetical 
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Figure 16. Graph for Showing Certain Programming 
Considerations 
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GH(joi) frequency response are assumed as A and B. The points a and g 
are the points where the margins occur. If they are perturbed so that 
the distances between them and the -1 4- jO point are increased, then 
they are allowed to move in any direction which has a positive dot 
product with vectors emanating from the -1 + jO point to these points. 
Suppose that a was perturbed in the direction 0 indicated in Figure 16. 
It is obvious that by moving a in this direction the vector from -1 + 
jO to a is increasing in magnitude. However, after a is perturbed it 
is no longer the point of interest. Some other point such as X is 
then the point under consideration, where X is in some neighborhood of 
a. From practical considerations it is known that if a moves in the 
direction 6 then a small neighborhood around a will move in the direc- 
tion 8. Let X be in this, neighborhood. The result is that X will be 
the new point of intersection with the real axis, and, furthermore, 
its distance from the -1 + jO point is less than what a's was. Similar 
results can be demonstrated for 6. 

These types of problems can be circumvented by perturbing a point 
with respect to two points in the complex plane. In fact consider the 
example in the last paragraph. Suppose that a is not only pushed with 
respect to the -1 + jO point, but it is also pulled with respect to 
the —7 — j 4 point. The permissible" region for the movement of a now 
becomes the intersection of the permissible region for pushing from 
the -1 + jO point and the permissible region for pulling with respect 
to the —7 — j 4 point. The result is the cross-hatched area in 
Figure 16. Movement of a anywhere in this region cannot result in 
the gain margin being decreased. 
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For the single channel example conditions did not exist to warrant 
pertubations with respect to more than one point. On the other hand 
the dual channel example required perturbing the aerodynamical gain 
margin with respect to two points. Runs in which this was not done 
resulted in a significant reduction in this margin. 

In neither example did the lead and lag phase margins of the first 
bending mode become active. In the single channel example, conditions 
just never prevailed. As for the dual channel example, conditions 
would have probably resulted if the magnitude of the first bending mode 
had not been controlled by the attenuation margin technique. Since the 
frequencies where these margins occur are very close to the frequency 
of the first bending mode, then it is quite natural that an increase 
in the first bending mode magnitude would have resulted in the reduc- 
tion of at least one of these margins. 

The program indicated for the dual channel example that better 
results could be obtained with one less zero in the numerator of the 
first channel's compensator and one less pole in the second channel's 
compensator. It did this by driving these to infinity. It also drove 
two poles in each channel to equal values. This probably indicates 
that if these poles were included in second order factors they would 
split into complex conjugates. However, the first order pole factors 
were chosen so that complex poles would not be allowed. 

One other fact which should be pointed out is that the program 
was used in the SIFR mode. However, because of the maximum step size 
choices (0.1 for the first 1000 iterations of the first example and 
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0.2 for the second 1000 Iterations and u.l for the second example) the 
program actually performed in the TIFR mode.* 

One phenomenon which should not pass without mention is the 
apparent unsmoothness of the convergence curves, Figures 7 and 15. In 
actuality, these curves should be discrete curves* For convenience 
they were drawn continuously. The sharp', abrupt changes, where the 
smallest stability margins make much greater gains than on other itera- 
tions, occur at iterations where the aerodynamical gain margin became 
inactive. This allowed the smallest stability margin to make a marked 
gain for one iteration. While this was occurring the aerodynamical 
gain margin was returning to activity. Once it became activeagain 
the rate of increase of the smallest stability margin decreased. On 
higher iterations the curve was smooth until the aerodynamical gain 
margin went inactive again, at which time the process was repeated. 

The overall effect of the program is a “ratchet" type, i. e. , once a 
margin is increased; it will not decrease. 


Of course again this is neglecting instances where constraints 
went from inactivity to activity. 
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CHAPTER VIII 


CONCLUSION, LIMITATIONS, AND SUGGESTED 
FUTURE STUDIES 

Summary 

In this dissertation, the theory for a compensator improvement 
algorithm has been presented. The goal from the onset was to accom- 
plish this by way of mathematical programming. Thus, in Chapter I 
a concise review of the more popular mathematical programming tech- 
niques was given. After this review a discussion of the uses of 
mathematical programming in the design of control systems was pre- 
sented. Also, a discussion of the uses of mathematical programming 
in the design of control systems was made. In this discussion it was 
pointed out that only a small amount of effort has been devoted to 
using mathematical programming as an aid in the design of control 
systems by classical means. Furthermore, it was shown that the tech- 
niques which have been developed suffer from some serious drawbacks. 
Thus, the thesis of this dissertation was to develop a computerized 
compensator design procedure* whictr circumvented these pitfalls. 

In Chapter II, some important- concepts for the measuring of 
expected performance of a control system were given. This involved 
defining relative stability in a way- somewhat different from the 
normal textbook definition. Also, concepts of relative attenuation 

and proper phasing were defined. Finally, using these the design of 
a compensator was formulated as a mathematical programming problem — 
which in the end resulted in a strict' constraint problem. 
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In Chapter III compensator' limitations for two possible iterative 
techniques for solving the problem formulated in Chapter II were pre- 
sented by the proving of two theorems. The first theorem showed that 
to be assured of being able to perturb n points in the GH(juj) plane 
in n optimal directions there must exist 2n coefficients for variance. 
On the other hand. Theorem 2 stated that if each point was given 180° 
of freedom for movement (a sub-optimal direction), then only n coef- 
ficients were needed for variance. From this it was deduced that a 
sub-optimal algorithm would be the most practical. 

. Then, in Chapter IV the development of a sub-optimal algorithm 
was made. ,The result was the evolvement of the constraint improve- 
ment algorithm. In this development" several definitions were given, 
e.g. , total improved frequency response, sum improved frequency 
response, improved solution, and active and passive constraint's. 

In order to employ the constraint improvement algorithm in 

Chapter IV, it was expedient to have the gradients of the active con- 

til 

Straints. These were found in Chapter V for a general j channel 
control system. Furthermore, the partials were derived so that push- 
ing .or pulling on points of the frequency response could be accom- 
plished with respect to any points desired in the complex GH(jto) 
plane. 

Next, the ideas and material in Chapters II, III, IV, and V were 
included in a computer program" called CIP (Compensator Improvement 
Program). In Chapter VI the general iterating procedure of this pro- 
gram was incorporated. In addition, several special programming 
techniques employed by CIP were presented in this chapter. 
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Chapter VII was used to demonstrate the practicality of CIP. This 
was illustrated by two large system' examples. These examples clearly 
showed the program's capability of handling single or multi-channel con- 
trol systems. A significant amount of improvement in the frequency 
response of both systems was seen after an application of CIP. 'Also, 
curves to show the convergence properties of CIP were given. In 
addition, several comments in regard' to proper specifications for the 
program were mentioned. 

Limitations and Concluding Remarks 

One of the limitations of CIP is that the initial compensator must 
be chosen to stabilize the system. This is the reason that the program 
was termed an "improvement program" rather than a design program. A 
major goal of the program is to improve stability margins, etc. , from 
one iteration to another. Obviously, if the system is initially un- 
stable then stability margins have no meaning. 

Another shortcoming of eiP is that a choice of the components of 
the c vector in Chapter IV must be'made. If the strict constraint 
problem has a solution which is reachable from the initial starting 
point, the choice of the c vector has little consequence other than to 
affect the rate of convergence. However, if the problem does not have 
an obtainable solution, then the choice of this vector will definitely 
determine the relative extremal where convergence occurs. Neverthe- 
less, it should be pointed out that if the initial guess at the 
solution is not a relative extremal then the solution at convergence 
will be better than the initial solution. 
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A very good property which CIF possesses is an inherent ability 
not to design an unstable' compensator, provided the step size is main- 
tained at a reasonable value. The reason for this is that CIP con- 
tinuously improves relative stability; thus the stability of the 
system cannot decrease. 

Although CIP requires a choice of the c vector - elements* it still 
has the capabilities of yielding a practical design on every run. As 
long as the input specifications of the program are properly made, 
the program cannot yield a compensator worse than the original compen- 
sator. CIP is not a design technique, but it is a design aid. 

Suggested Future Studies 

There are several areas in which the work in this dissertation 
can be extended. One such study could involve using the constraint 
improvement algorithm In other design problems Ln engineering and 
science. This author does not see any reason that it could not be 
used to make improvements in any design 'where the number of variables 
is greater than the numbei of constraints to be controlled and where 
the gradient vectors of the constraints are deterministic. 

Also, it is foreseen by this author that the constraint improve- 
ment algorithm could be the basis of a new or extended gradient algo- 
rithm for nonlinear programming. For example, if any of the elements 
of the c vector are set to zero then the determined directional vector 
will lie in the - tangent planes of the constraints corresponding to the 
c's with zero value. Of course this would be similar to the gradient 
projection technique mentioned 'in Chapter I; However, it is deemed by 
this author that by using the constraint improvement approach an 
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optimal gradient projection algorithm can be developed. Up to the 
present such an algorithm has not been developed. 

In regard to future studies in compensator design, there would be 
nothing wrong with starting with physical electrical networks, rather 
than transfer functions. If a program started with a network and varied 
the elements for making the improvements described previously, the end 
results would be the actual network needed. The practicality of this 
network would depend upon the constraints placed on the network ele- 
ments . 

A compensator design procedure could be devised using the con- 
straint improvement algorithm on the Routh-Hurwitz array. By forming 
the characteristic equation as a function of the compensator coeffic- 
ients, the first two rows of the Routh array can be formulated as 
functions of these compensator coefficients. Since it is known how the 
other rows of the array are formed from the first two rows, the changes 
in the elements of the first column of the array with respect to the 
compensator coefficients could be determined by an application of the 
chain rule for partial derivatives. Then, the constraint improvement 
algorithm could be used to drive all the negative elements of the first 
column positive, as long as the number of negative elements did not 
exceed the number of compensator coefficients. If all the elements are 
driven positive then a certain amount of relative stability could be 
achieved by evaluating the characteristic equation at (s + a) where a is 
a positive real number; the previously mentioned procedure can now be 
applied to the new characteristic equation. If in this application all 
elements of the first column could again be driven to positive values, 
then it would be known that no pole of the closed loop system has a 
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real part greater than - a. This process coaid be repeated until a 
desired value of a is achieved or until all the elements of the first 
column of one of the characteristic equations cannot be driven 
positive. 
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APPENDIX A 


COMPENSATOR IMPROVEMENT PROGRAM 

In the following is a complete Fortran version of the Compensator 
Improvement Program. The program is completely self-contained, i.e., 
it does not require any system library, etc. The necessary input to 
the program is explained in the comment statements at the beginning of 
the main program. Furthermore, all inputs except the frequency 
response points are printed out with explanations of the input speci- 
fications, The other output is, also, explained by certain comments 
printed out with the information. 
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oooooooooooooooooooooooociooooooooonociooonoooooooon 


MAIN PROGRAM 

definitions OF I/O variables 

KCHNL -NO. OF CHANNELS FED BACK 

NUMC(I) -NO. OF COMPENSATORS IN I-TH CHANNEL 

NRATOR(I,J) -numerator order of j-th compensator in the i-th channel 
Ndenom ( I » J) -denominator order of j-th comp, in i-th channel 

XCOMN ( I » J) -NUMERATOR COEFFICIENTS OF J-TH COMP. IN I-TH CHNL. 
YCOMN(l.J) -DENOM. COEFFICIENTS OF J-TH COMP. IN I-TH CHNL. 

OMEGA { i ) -I-TH FREQ. (ASSUMED TO BE IN HZ.) 

GRA(I.J) -J-TH REAL PART OF OPEN LOOP FREQ. RESP. OF I-TH CHNL. 

GIA(Irj) -J-TH IMAG. PART OF OPEN LOOP FREQ. RESP. OF I-TH CHNL. 

KSTART -STARTING ITERATION NO. 

KQUIT -STOPPING ITERATION NO. 

KPOINT -NO. OF POINTS FrOM OPEN LOOP FREQ. RESPONSE USED 
KPRINT - NO. OF ITERATIONS SKIPPED BETWEEN PRINTING OF INFOR. 
STFMAX -MAXIMUM CHANGE TO BE MADE IN COMPENSATOR COEFFICIENTS 
ON ANY ONE ITERATION (PROBABLY NO MORE THAN 30* OF THE 
SMALLEST COMPENSATOR COEFFICIENT OF THE INITIAL 
COMPENSATOR) 

STPMIN - MINIMUM STEP SIZE DESIGNATOR 
Flo & Fll - FREQUENCIES BETWEEN WHICH G.M.'S ARE FOUND 

F12 & F13 - FREQUENCIES BETWEEN WHICH P.M.'S ARE FOUND 

FMIN - A.M.'S ARE FOUND FOR FREQS. ABOVE THIS FREQ. 

variables for gain margin radii designations 

IF FREQ. .LE. FI DESIRED MARGIN = R1 

IF FREQ. .GT. FI BUT .LT. F2 DESIRED MARGIN = R2 

IF FREQ. .GE. F2 DESIRED MARGIN = R3 

VARIABLES FOR PHASE MARGIN RADII DESIGNATIONS 
IF FREQ. .LE. F3 DESIRED MARGIN = R4 

IF FREQ. .GT. F3 BUT .LT. F4 DESIRED MARGIN = R5 

IF FREQ. .GE. F4 DESIRED MARGIN = R6 

VARIABLES FOR STABILITY MARGIN RADII DESIGNATIONS 

IF FREQ. .LE. F5 DESIRED MARGIN = R7 

IF FREQ. .GT. F5 BUT .LT. F6 DESIRED MARGIN = R8 

IF FREQ. .GE. F6 DESIRED MARGIN = R9 

VARIABLES For ATTENUATION margin radii designations 
IF FREQ. .LE. F7 DESIRED MARGIN = RIO 

IF FREQ. .GT. F7 BUT .LT. F8 DESIRED MARGIN = Rll 

IF FREQ. .GE. F8 DESIRED MARGIN = R12 

GAINU)-DEN0TES INITIAL D. C. GAIN VALUE FOR I-TH CHANNEL 

KNR(I) -NUMBER OF NUMERATOR. COEFS. FOR I-TH CHANNEL 
KDR(I) -NUMBER OF DENOM. COEFS. IN I-TH CHANNEL 
KONT ( I ) -D.C , DESIGNATOR FOR I-TH CHANNEL 

KONT(I)=l GAIN ALLOWED TO VARY 
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KONT ( I ) =2 GAIN NOT ALLOWED TO VARY 

Klf-M -NO. CHANNELS THAT FREQ. RESP. INFIRmATION IS TO BE READ IN 

PPT(I) -POINTS THAT THE CRITICAL FREQUENCIES WILL BE 
PERTURBED WITH RESPECT TO (COMPLEX POINTS) 

1=1 GAIN MARGIN POINT 
1=2 PHASE MARGIN POINT 
1=3 stability MARGIN POINT 
1=4 ATTENUATION MARGIN POINT 

lsn(i) - denotes whether points are to be pushed or pulled 

LSN=-1 POINT TO BE PULLED 
LSN=+1 POINT TO BE PUSHED 

INCGMS -INDICATES WHETHER G.M.*S ARE TO BE ARTIFICALLY 
INCLUDED AS S.M.’S 
INCGMS=0 NOT INCLUDED 

INCGMS=1 INCLUDED 

INCPMS -INDICATES WHETHER P.M.»S ARE TO BE ARTIFICALLY 
INCLUDEED AS S.M.»S 
INCPMS=0 NO INCLUDED 

INCPMS=1 INCLUDED 

some Interior variable definitions 


GCR(IrU) 

GC I ( I . u ) 
GCOMR ( 1 » U) 
GCOMR ( I # J) 
GR(I) 

GI(I) 


-REAL PART OF COMPENSATOR FREQ. RESP. AT SOME ITERATION 
-IMAG. PART OF COMPENSATOR FREQ. RESP. AT SOME ITERATION 
-REAL PARTS OF I-TH CHNL* OPEN LOOP FREQ. RESp. 

-IMAG. PARTS OF I-TH CHnL. OPEN LOOP FREQ. RESP. 

-REAL PARTS OF TOTAL OPEN LOOP FREQ. RESP. 

-IMAG. PARTS OF TOTAL OPEN LOOP FREQ. RESP. 


***** THERE ARE 13 READ STATEMENTS ***** 

DIMENSION XCOMNU0.50) »YCOMN(10»50> rPRY(50> #PRX(50) .STBM(99) » 

1 PX ( 50 > » PY ( 50 ) . RQ ( 99 ) . GR < 999 ) . GI ( 999 ) » OMEGA ( 999 ) . GRA < 5 . 999 ) » 

2 GIA(5*999) »G(20»99) »DV(50) ► WEIGHT (50) »BCOMN( 10 .50 ) , 

3 6C0MD ( 10.50 ) »GCR(5» 999) »GCI (5»999) .GCOMR (5. 999) . 

4 GCOMI (5.999) . NUMC ( 20 ) .NRATOR ( 10 . 20 ) . NDENOM( 10 . 20 ) . CNUM( 10 ) » 

5 CDOM(IO) .KNR(IO) »KDR(10) »COTN(10.50) »COTD( 10.50) 

DIMENSION KACT(99) »SML(99) 

DOUBLE PRECISION G.DV. WEIGHT 
DIMENSION KONT (20 ) » KPTS(99)» GAlNdO) 

DIMENSION TYPE(99> 

DIMENSION PPT (4 ) » LSN(4) 

COMMON type 

integer type 

Complex PPT 

READ (5» 5) KCHNL 

PEaD ( 5» 5) (KONT(I).I = l»KCHNL) 

READ(5»5) (NUMC(I) »I=1»KCHNL) 

WRITE (6.1) KCHNL 

i Format < *o* »sx . ’number of channels fedback=*,i5) 

WRITE(6.3) ( KONT (I) .1=1. KCHNL) 

3 FORMAT('0*»5X»»D.C. GAIN CONSTRAINT DESIGNATOR FOR EACH CHANNEL ( 
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1K0NT=1* ALLOWED TO VARY! K0NT=2. HELD CONSTANT ) * /6X*8 ( 12* 10X) ) 
WRITE (6*4) CNUMCCI) . I=1*KCHNL) 

• .4 FORMATC *0» *5X* 'COMPENSATORS PER CHANNEL* *1015) 

DO 2 I = 1*KCHNL i . 

KNAT=NUMC(I) 

REaDC5*5> (NRATOR(I*J) *J=1*KNAT) 

WRITEC6.6) I» (NRATOR(I'J) »J=1»KNAT) 

. , ; 6 FORMAT ('O'* 5X'»» CHANNEL NO. » *I2*2X* 'NUMERATOR ORDERS' *2X. 1015) 

; , • : READ (5* 5) (NDENOM<I*J)'J=l*KNAT) 

WRITE(6*7) I* (NDENOM(I*J) *J=1*KNAT> 

7 FORMATC « *5X* 'CHANNEL NO. ♦ f I2*2Xr 'DENOMINATOR ORDERS' * 1015) 

2 CONTINUE 

, .5 FORMAT ( 1615). * . 

• ~.v REAO(5*10> KSTART*KQUIT*KIFMfKPOlNT*KPRINT» R1 *F1*R2*F2*R3* 

1 R4*F3*R5fF4.R6* R7* F5*R8*F6*R9* RlOf F7f Rllf F8*R12* FMIN*F10* 

2 Fllf F12»F13f STPMAXf STPMIN 

10 FORMAT (5I5/5F10 .5/5F10.5/5F10 .5/5F1 0*5/ 6Fl 0.5) 

WRITE (6 * 11 )KSTART *KQUIT f KIFMf KPOINT * KPRINT 

11 FORMAT l 'O'* IX* 'START ITeR.=* r I5*2X* 'STOP ITER»=* * I5*2X* 'No. CHNL. 
1FREQ. RESP. IN=* *I5*2X*'N0. OF FREQ. POINTS=* *I5»2*' 'PRINT INCREME 
2NTs**I5) 

WRITE(6*25)STPMAX*STPMIN 

25 FORMAT( 'O' *5X* 'MAXIMUM DESIGNATED STEP SIZE = » *F10.5/6X* 'MINIMUM D 
1ESIGNATED STEP SIZE ='*F10.5) 

WRITE (6 * 12) 

12 F0RMAT('0»5X* 'DESIRED GAIN MARGIN RADII DESIGNATIONS’) 

WRITE (6* 13) F1.R1* F1*F2*R2* F2*R3 

13 FORMATC *0»»5X* 'IF FREQUENCY .LE. * *F10.5*5X* 'DESIRED MARGIN IS'* 

1 F10.5/6X* ' IF FREQUENCY .GT. ' *F10.5*2X* 'BUT .LT. » *Fl0.5*2X* 'DESIRE 
2D MARGIN IS* *F10.5/6X* » IF FREQUENCY .GT. » *F10.5*2X* 'DESIRED MARGI 
3N IS' *F10.5) 

WRITE (6* 17) F10 *F11 

17 FORMATC »*5X*'GAIN MARGINS ARE DETERMINED BETWEEN THE FREQUENCIES 
1 OF* *F10«5*2X* ' AND* * F10 .5) 

WRITE(6*14) 

14 FORMATC 'O' *5X* 'DESIRED PHASE MARGIN RADII DESIGNATIONS') 

WR1TEC6* 13) F3*R4# F3*F4*R3* F4*R6 

WRITEC6* 18) F12*F13 

18 FORMATC* »*5X* 'PHASE MARGINS ARE DETERMINED BETWEEN THE FREQUENCXE 
• IS OF»*F10.5*2X*'AND'*F10.5> 

WRITE (6* 15) 

15 FORMATC »0'*5X*'DESIRE0 STABILITY MARGIN RADII DESIGNATIONS') 

WRITE (6* 13) F5*R7* F5*F6*R8* F6*R9 •' 

WRITE(6*16> 

16 FORMATC »0"5X* 'DESIRED ATTENUATION MARGIN RADII DESIGNATIONS') 
WRITEC6*13) F7*R10* F7»F8*Rll* F8*R12 

WRITE (6* 19) FMIN 

19 FORMATC '*5X*' ATTENUATION MARGINS ARE FOUND FOR FREQS* ABOVE' » 

1 F10.S) 

READ (S' 50) CGAINCD flsliKCHNL) 

WRITE(6*20) ( X * GAIN ( I ) *I=1*KCHNL) 

20 FORMATC '0"5X *2 ( 'CHANNEL NO. ’* I3*1X* ' INITIAL D.C. GAIN IS'*F10.5* 
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1 5xn 

READ<5*50) (PPT(I) »I=1»4) 

WRITE(6#22) (PPT(I) *1=1*4) 

22 format ( ’ o • # sx # » pertubation points for gain# phase# stability# and 

1ATTENUATION MARGINS# RESPECTIVELY: */6X#4(*ReAL* #F6.2»2X» 'IMAG. * # 

2 F6.2#2X>) 

READ (5#5) (LSN(I) »I=l»4) 

WRITE(6»23) (LSN(I) #I=1»4) 

23 FORMAT! *0’#5X#’0EN0TING WHETHER EACH OF THE PRECEDING POINTS ARE * 

1 #'TO BE PUSHING OR PULLING PoINTS(PUSHING=+1 # PULLING=-1)* /6X» 

2 4(12, 10X)) 

READ(5#5) INCGMS# INCPMS 
WRITE (6#24) INCGMS# INCPMS 

24 FORMAT( ’O’ #5X# ’DENOTING WHETHER GAIN OR PHASE MARGINS ARE ARTIFICA 

ILLY INCLUDED AS STABILITY MARGINS (NOT INCLUDED=0# INCLUDED=1 ) ’/6X# 
2 2(I2#10X)) 

KVARY=0 

DO 21 K=1#KCHNL 
LAMP=NUMC(K) 

KNR(K)=0 
KDR(K)=0 
DO 21 I=1#LAMP 
KVARY=KVARY+NRATOR (K # I ) 

KVARY=KVARY+NOENOM (K # I ) 

KNR(K)= KNR(K) + NRATOR(K#I) + 1 
21 KDR(K)= KDR(K) + NDENOM(K#I> + 1 
DO 29 I=1#KCHNL 

29 IF(KONT(I) .EQ.1>KVARY=KVARY+1 
DO 42 K=1#KCHNL 
LNC= KNR(K) 

LDC= KDR(K) 

READ (5 #50) (XCOMN (K # I ) # I— 1 # LNC ) 

42 REAO ( 5# 50 ) ( YCOMN(K # I ) # 1 = 1 # LDC ) 

50 FORMAT (8F1Q • 5 ) 

60 FORMAT('0*#6X#»INITIAL COMPENSATOR COEFFICIENTS’) 

WRITE(6#60) 

DO 72 K=1#KCHNL 
LNC=KNR(K) 

LDC= KDR(K) 

WRITE(6#62)K 

62 FORMAT( *0* »5X» ’CHANNEL NO. * #12# 2X#' COMPENSATORS - FACTORED FORM’) 
WRITE (6 #68) 

68 FORMATt ’O’ #5X# ’NUMERATOR COEFFICIENTS’) 

WR1TE(6#70) ( XCOMN (K# I ) # 1=1# LNC > 

WRITE (6#69) 

69 FORMAT! ’0’#5X# ’DENOMINATORS COEFFICIENTS*) 

72 WRITE(6#70) (YCOMN(K#I> # I=1#LDC) 

70 FORMAT ( ’ • #5X#10F10.5) 

C MODIFICATION OF FREQ. RESP. INFOR. BY CONTANT COMPENSATOR 
00 135 J=1»KIFM 

135 READ(5#140) (OMEGA ( I ) #GRA( J# I ) #GIA( J# I ) # 1-1 #KPOINT) 

140 FORMAT (9F8. 5) 


91 



IF(KIFM.GE.KCHNL)GO TO 150 

K= KIFM + 1 

DO 148 J=K»KCHNL 

DO 148 I=l»KPOINT 

XV= OMEGA ( I ) * 6.2831853 

GRA(J#I)= -OMEGAU) * XV * GIA(J-Ifl) 

148 GIA(J»I)= OMEGA(I) * XV * GRA(J-lrl) 

150 CONTINUE 

DO 149 J=1»KCHNL 

DO 149 1 = 1 • KPOINT 

GRACJ»I)= GRA(J.I) * GAIN ( J) 

149 6IA ( J» I )= GIA ( J» I ) * GAIN(J) 

190 CONTINUE 

DATA STEP # KHOP t SML2 1 PSQL * SBC2/1 . 0E-02 1 0 1 0 . 0 r 1 . 0E+20 >0.0/ 
11=0 
12=0 

DO 195 K=1»KCHNL 
11= KNR(K) ♦ II 
195 12= KDR(K) + 12 
LOX= KSTART 
200 CONTINUE 

LPR£SV=KVARY 

NM=0 

C EVALUATION OF VARIABLE COMPENSATOR AT CHOSEN FREQS. 

DO 210 K=1 » KPOINT 
GR (K ) =0 . 0 
GI (K)=0,0 

XV=OMEGA<K) *6.2831853 
DO 209 1=1 »KCHNL 
KCOtMP= NUMC(I) 

lnot=o 

knot=o 

GCR(I»K)= 1.0 
GCl ( I #K)= 0.0 
DO 208 J=1 » KCOMP 
NTR= NRATOR(I»J)+l 
NTU= NDENOM(I» J>+1 
DO 204 M=1 »NTR 

204 CNUM(M)= XCOMN(I»M+KNOT) 

DO 205 M=1»NTD 

205 CDOM(M)= YCOMN ( I » M+LNOT ) 

KNOT= KNOT + NTR 

LNOT= LNOT + NTD 
K2= NTR-1 
K3= NTD-1 

CALL P0LFV(CNUM.K2»XV.CnR»CNI) 

CALL P0LFV(CD0M.K3»XV»CDRrCDI) 

CD= CDR**2 ♦ CDI**2 
ACR=GCR(I»K> 

ACI= GCI(I»K) 

ACOMR= (CNR * CDR + CNI * CDD/CD 
ACOMI=(-CNR * CDI + CNI * CDRJ/CD 
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GCR(I»K)= ACR * ACOMR - ACI * ACOMI 

208 GCI<I#K)= ACR * ACOMI + ACI *■ ACOMR 
GCOMR(I»K)= 6RA(I»K)*GCr(I»K) - GIA(I»K)*6CI(I»K) 
6COMI ( I »K) = GRA ( I # K ) *GC I ( I » K ) + GIA( I » K ) *GCR ( I » K ) 
GR(K)= GR(K) + GCOMR(I»K) 

209 GI(K)= GI(K) + GCOMIU’K) 

210 CONTINUE 

C DETERMINATION OF GAIN MARGINS POINTS BETWEEN F1‘ AND F2 

call g ainmg ( gr » g i » kpoi nt » nm » fi o » f i i , kpts » stbm Iomega ) 
ngms=nm 

c setting desired stability radii of g.m.*s 

KPM=NM+1 

IF (NM.EQ » 0 ) GO TO 213 
DO 212 I=1»NM 
i TYPE(I)= *G* 

KWHICH=KPTS(I) 

FREHZ=OMEGA ( K WH I CH ) 

IF(FREHZ.LE.F1)RQ(I)=R1 

IF(FREHZ.GT,F1)RQ(I)=R2 

IF(FREHZ.GE.F2)RQ(I)=R3 

212 CONTINUE 

213 CONTINUE 

C DETERMINATION OF P.M. BETWEEN F3 AND F4 

CALL PHASEM ( GR » G I f KPO INt » NM r F12 r F 1 3 » KPTS r STbM r OMEGA ) 
IF(NM.LT«KPM)GO TO 215 

c setting desired stability radii of p.m.*s 

DO 214 I=KPMfNM 
TYPE(I>= »P* 

KWHICH=KPTS(I) 

FREHZ=OMEGA (KWHICH) 

IF(FREHZ.LE.F3)RQ(I)=R4 

IF(FREHZ.GT.F3)RQ(I)=R5 

IF(FREHZ.GE.F4)RQ(I)=R6 

214 CONTINUE 
KPM=NM+1 

215 CONTINUE 
IF(NM.EQ.O)GO TO 221 

klast=nm 

DO 220 I = 1 » KLAST 

IF( (I.LE.NGMS) .AND. (INCGMS.EQ.l) >GO TO 2l9 
IF( (I.GT.NGMS) .AND. ( INCPMS.EQ.1) )GO TO 2l9 
GO TO 220 

219 KPM=KPM+1 

nm=nm+i 

KPTS(NM)=KPTS(I) 

STBM t NM ) =ST BM ( I ) 

RQ(NM)=RQ(I) 

TYPE(NM)=*S« 

220 CONTINUE 

221 CONTINUE 
KSTBM=KPM 
RPT=1.0 


93 



NS6=1 

FQMIN=0.0 

C DETERMINATION OF STABILITY MARGINS . • . 

CALL SRMlNS(GR»GI»KP0INT»NM»RPT»NS6»FQMIN»KpTS»STBM*0MEGA) 

C SETTING DESIRED STABILITY MARGINS 
IF<NM.LT.KPM)GO TO 216 
DO 230 I=KPM»NM 
TYPE ( I ) = • S* 

KWHICH=KPTS(I) ? 

FR£HZ=OMEGA ( KWH I CH ) 

IF(FREHZ.LE.F5)RQ(I)=R7 
IF(FREHZ.GT.F5)RQ(I)=R8 
IF(FREHZ.GE.F6)RQ(I)=R9 
230 CONTINUE 

C CHECKING TO SEE IF ANY P.M.«S. G.M.*S» OR S.M.’S ARE EQUAL 
C IF THERE RESULTS SOME THAT ARE EQUAL ONLY THE FIRST IS RETAINED. 

DO 228 LB=2,KSTBM ' : 

DO 228 I=KSTBM#NM 

IFlKPTS(LB-l) .NE.KPTSd) )GO TO 228 

NM=NM-1 

DO 226 L=I»NM 

KPTS(L)= KPTS<L+1) ' 

STfciM (L) = STBMU.+1) 

RQIL)= RQIL+l) 

226 TYPE (L) = TYPE(L+1) 

228 CONTINUE 
KPM=NM+1 

216 CONTINUE 

KMIN=NM , . 'r 

RPT=0.0 ' ■ 

NSG=-1 • 

FqmIN=FmIN 

C DETERMINATION OF ATTENUATION MARGINS . . 

CALL SRMINS (GR » G I » KPO INT » NM » RPT . NSG * FQM IN » KpTS » STBM , OMEGA ) 

C SETTING DESIRED ATTEN. MARGINS 

IF(NM.LT.KPM)GO TO 217 . > ! 

DO 232 I=KPM»NM 
TYPE<I)= • A* 

KWHICH=KPTS(I) 

FREHZ=OMEGA (K WHICH) 

IF(FREHZ.LE.F7)RQ(I)=R10 
IF IFREHZ.GT .F7)RQ(I)=Rll 
IF(FREHZ.GE.F8)RQ(I)=R12 
232 CONTINUE ; 

217 CONTINUE ( 

SBC1=R1 - 

C DETERMINING SMALLEST STABILITY MARGINS OF PRESENT ITER. AND ALL ITER. 
SML1= 100.0 
DO 290 I=1»KMIN 

IF (STBM ( I ) .GT.SML1) GO To 288 1 r 

SML1= STBM ( I ) 

288 CONTINUE 
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IF(STBM(I) .GT.SBCD60 TO 290 
SBC1= STBM(I) 

290 CONTINUE 

IF ( SBC2 • GE t SBC1 ) GO TO 298 

SBC2= SBC1 

IBEST= LOX 

DO 292 K=l* KCHNL 

LNC= KNR(K) 

LOC= KDR(K) 

DO 291 I=1»LNC 

291 BCOMN(K*I)= XCOMN(K#I) 

DO 292 I=1*LDC 

292 BCOMD(K»I)= YCOMN(K»I) . 

.298 CONTINUE 

c checking satisfaction OF SYSTEM REQUIREMENTS 
DO 320 1=1 rNM 
PORM= 1.0 

IF l I ,GT.KMlN)PORM=-l . 0 
310 IF { (STBM( I )-RQ ( I ) ) *PORM) 350*320 r 320 
320 CONTINUE 

WRITE(6»330) 

330 format ( • o * * isx » * ***** all system requirements have been met ****** 
l) ' 

340 CONTINUE 

CALL OTPT1 ( STBM . OMEGA * KpTS * NM t XCOMN * YCOMN r KM IN » RQ * LOX , KCHNL * NUMC * 

1 NRAT0R»NDEN0M#PRX*PRY»I1»I2) 

WRITE(6i341) IBEST 

341 FORMAT(«0*»5X» ****** BEST COMPENSATORS WITH RESPECT To STABILITY * 
l****t//6Xr ’OCCURRED ON STEP* , I4»2X* ’AND THEIR COEFFICIENTS ARE ! * } 

CALL MULOUT { KCHNL * NUMC * NR ATOR * NDENOM * KNR * KDR * BCOMN * BCOMD > 
WRITE<6*345) SBC2 

345 F0RMAT(*0»»21Xr ‘SMALLEST STABILITY MARGIN FOR THE BEST COMPENSATOR 
/ =*»F10.8) 

WRITE(6»347) 

347 FORMAT (»0*»5X» ****** COMPENSATORS AND COMPENSATED FREQUENCY RESPON 
1SE AT THE LAST ITERATION PERFORMED ARE AS FOLLOWS ******) 

CALL MULOUT { KCHNL » NUMC * NR ATOR r NDENOM t KNR t KDR » XCOMN * YCOMN ) 
WRITE(6.346) 

346 FORMAT(*0**9X» 'COMPENSATED FREQUENCY RESPONSE*//10X» 'FREQUENCY* * 

1 2X *' MAGNITUDE '»3X»' ANGLE') 

DO 349 I=l»KPOINT 

GMTE= SQRT(GR(I)**2 + GI(I)**2) 1 

AGLE= ATAN2(GI(I) »GR(I) )*57.3 
WRITE (6. 348) OMEGA ( I ) » GMTE * AGLE 

348 FORMATt* • *7X»F10.5»lXrF10.5rlXfFl0.5) 

• 349 CONTINUE 

STOP 

350 CONTINUE 

C STEP SIZE SELECTING 

IF(L0X.EQ.KQUIT)WRITE(6»351) 

351 F0RMAT('0*»5X» ****** TERMINATION REASON - MAXIMUM ITERATIONS ***** 
1*) 
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IF(LOX.EQ.KQUIT)WRITE<6*400)STEP 

IFILOX.EQ.KQUITJGO to 340 

IF(LOX.EQ.KSTART)GO TO 354 

ADD=0.0 

MAD=0 

P0RM=1.0 

DO 355 1=1 »NM 

IF(I.GT.KMlN)PORM=-1.0 

IF(P0RM*(STBM(I)-RG(I) ) .GE.O.O)GO TO 355 
DO 352 J=1#NML 

352 IF(KPTS(I).£Q.KACT(J))GO TO 353 
MAO=MAD+l 

GO TO 355 

353 CONTINUE 

C IF IT IS DESIRED TO HAVt ALL CONSTRAINTS TO BE IMPROVED AT EVERT 
C ITERATION REMOVE THE C FROM COLUMN 1 OF THE FOLLOWING CARD 
C IFlPORM*(STBMlI)-SML<I> > .LT .-1 . OE-O5>60 TO 360 
ADO=AOD+PORM*(STBM(I)-SML(J) ) 

355 CONTINUE 

IF(MAD.EQ.NML)ADD=1.0 
IFtADD.LE. 0.0)60 TO 360 

354 CONTINUE 
GO TO 371 

360 STEP= STEP/2.0 

IF(STEP.LT.STPMIN )WRITE(6*365>STPMIN 
365 FORMAT! *0'*5X #'****♦ TERMINATION REASON - STEP SIZE IS LESS THAN * 
1 *F10>5’2X> •*****! ) 

IFtSTEP.LT.STPMIN )G0 TO 340 
LGX= LOX - 1 
GO TO 450 

371 STEP=1. 41416 * STEP 
373 CONTINUE 
SML2=SML1 

. IFtSTEP.GT.STPMAX)STEP= STPMAX 
C OUPUT CONTROL 

IF (KHOP«GT • 1 ) GO TO 410 
KH0P=KPRINT 
WR1TE16.400) STEP 

400 format ( * o* * isx» 'present step size =*»fio.7> 

CALL OTpTl ( STBM t OMEGA * KpTS » NM t XCOMN t YCOMN » KMIN r RQ » LOX r KCHnL t NUMC • 

1 NRAT0R»NDEN0M»PRX»PRY»I1»I2) 

60 TO 420 
410 KH0P= KHOP - 1 
420 CONTINUE 

C SELECTING ACTIVE CONSTRAINTS 
K=0 

DO 411 1=1 * NM 

IF(I-1.EQ*KMIN>KMIn=K 

P0RM=1.0 

IF(I.GT.KMlN)PORM=-1.0 
IF(PORM*STBM(I).GT.PORM*RQ(I))GO TO 411 
K=K+1 
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KPTS<K>= KPTSU) 

TYPE(K)= TYPE ( I ) 

SML(K>= STBmCII 
KACT(K)= KPTS(I) 

411 CONTINUE 
NM=K 
NML=NM 

C CALCULATION OF GRADIENTS OF ACTIVE CONSTRAINTS 
RPT=1.0 

CALL PARCLT ( XCOMN * YCOMN , GR » GI * OMEGA t NM t NRATOR • NDENOM • 

1 KCHNL * NUMC * KONT i GCOMR » GCOMI »G*PPT • LSN * KPARC * KPTS * KnR * KQR ) 

C SET OOT PRODUCT VECTOR 
DO 422 K=1*NM 

422 WEIGHT(K)=1.0 

C CALCULUTE DIRVECTIONAL VECTOR 
LR£=0 
KR£=o 

423 IFINM.GT.LPRESV)WR1TE<6.415) 

IF(NM.6T.LPRESV)G0 TO 340 

CALL OIRVEC (G*NM*KP ARC*DV* WEIGHT) 

415 FORMAT! »0»*5X* ****** TERMINATION REASON - NO* OF ACTIVE CONSTRAINT 
IS IS GREATER THAN THE No. OF ALLOWABLE VARIABLES *«■***’) 

DO 426 1=1*11 

426 PRX ( I ) = DV(I) 

DO 427 1=1*12 

427 PRY ( I ) = DV ( Il+I ) 

1F(KRE.EQ.1)G0 TO 433 

C CKECKING POSSIBLE NEGATIVENESS OF ANY COMPENSATOR COEF. 
IF(LRE.GE.I1+I2)G0 TO 433 
LRE=LRE+1 
K2=0 
K3=0 

DO 431 K=1*KCHNL 
LNC=KNR(K) 

LDC=KDR(K) 

DO 429 I=1*LNC 
K2=K2+1 

IFlXCOMNiK*I) «GT.1.0E-05)G0 TO 429 
IF(PRX(K2) .GE.O.OJGO TO 429 

lpresv=lpresv-i 

kre=i 

DO 428 J=1 *NM 

428 G(J,K2)=0.0 

429 CONTINUE 

DO 431 I=1»LDC 
K3=K3+1 

IF(YCOMN<K*I) .GT.1.0E-05>GO TO 431 
IF(PRY(K3) .GE.O.O)GO TO 431 

lpkesv=lpresv-i 

kre=i 

DO 430 J=1*NM 

430 G(j,H+K3)=0.0 
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431 CONTINUE 

IFCKRE.EQ.DGO TO 423 
433 CONTINUE 
PSQ= 0.0 
DO 438 I=1»I1 
PX(I)= PRXM) 

438 PSQ= PSQ + PX ( I ) **2 
DO 440 I = 1 t 12 
PY(I)=PRY(I) 

440 PSQ= PSQ + PY ( I ) **2 
PMG= SORT (PSQ) 

PSQL= PSQ 
DEL= STEP/PMG 
DO 462 K=1»KCHNL 
LNC= KNR(K) 

DO 462 I=1>LNC 
462 COTN(K#I)= XCOMN<K»I) 

DO 464 K=1»KCHNL 
LDC= KDR(K) 

DO 464 I=1»LDC 

464 COTD(K»I>= YCOMN(K» I ) 

GO TO 465 

450 DEL= DEL/2.0 

DO 467 K=1»KCHNL 

LNC= KNR(K) - ' 

DO 467 I=1»LNC 

467 XCOMN ( K » I ) — COTN<K»I> 

DO 468 K=1*KCHNL 
LDC= KDR(K) 

DO 468 I=1»LDC 

468 YCOMN(K,I)= COTD(K.I) 

465 CONTINUE 
KKK=0 

DO 470 K=1*KCHNL 
LNC= KNR(K) 

DO 470 I=1»LNC 
KKK= KKK+1 

XCOMN (K . I ) — XCOMN (K . I ) + DEL * PX(KKK) 
470 IF(XCOMk(K»I) .LT.O.O)XCOMN(K»I)=0.0 
KKK=0 

00 490 K=1»KCHNL 
LDC= KDR(K) 

DO 490 I=1*LDC 
KKK= KKK+1 

YCOMN(K»I)= YCOMN (K » I ) + DEL * PY(KKK) 
490 IF ( YCOMN (K* I ) *LT • 0 • 0 ) YCoMN (K r I ) =0 .0 
500 CONTINUE 

LOX= LOX + 1 
60 TO 200 
END 
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SUBROUTINE PARCLT ( XCOMN , YCOMN » GOR » GOI » OMEGA , NFREQ * NRATOR » NDENOM , 

1 KCHNLfNUMC»KONT»GCOMR»GCOMI»G»PPT»LSN»NPARC»KPTS»KNR»KDR) 

C 

C PROGRAM FOR CALCULATING THE CHANGE of A FREQUENCY RESPONSE WITH 
C RESPECT TO A COnPENSATOR COEFFICIENTS 
C 

c definitions OF I/O VARIABLES 
c 

C XCOMN { I »U) -NUMERATOR COEFs. OF COMPENSATOR IN I-TH CHANNEL 
C YCOMN(I,J) -DENOM. COEFS. OF COMPENSATOR IN I-TH CHANNEL 
C GOR Cl) -I-TH REAL P ART OF OPEN LOOP FREQ. RESP. 

C GOI (I) -I-TH IMAG. PART OF OPEN LOOP FREQ. RESP. 

C OMEGA ( 1 ) -I-TH FREQUENCY RESPONSE POINT 
C NFREQ -NUMBER OF MARGINS TO BE IMPROVED 

C NRATOR ( I » J) -NUM. ORDER OF J-TH COMP. IN I-TH CHANNEL 
C NDENOM ( I » J) -DEN. ORDER OF J-TH COMP. IN I-TH CHANNEL 
C KCHNL -NUM. OF CHANNELS 

C NUMC(I) -NUM. OF COMPS. IN I-TH CHANNEL 
C KONT(I) -GAIN CONTROL NUM. FOR I-TH CHANNEL 

C GCOMR( I » J) -REAL PART OF J-TH CHANNEL COMP. FREQ. RESP. AT J-TH FREQ. 
C GCOMI(IfJ) -IMAG. PART OF J-TH CHNL» COMP. FREQ. RESP. AT J-TH FREQ. 

C GU»J) -J-TH PARTIAL OF I-TH FREQ. 

C Z -NEG. OF POINT FOR WHICH PARTIALS ARE DESIRED 

C L -NO. OF POINTS TO TREAT AS STABILITY MARGINS<THE REMAINING 

C ARE CONSIDERED AS ATTENUATION MARGINS) 

DIMENSION C(10)tD(10)'£(10)rGR(50)rGl( 50 )» OMEGA ( 999) r Y ( 10 ) t X (10 ) r 

1 NUMC(20) rKONT(lO) »G(20»99) *G0R(999) rGOI (999) *NRATOR(10»20 ) • 

2 NDENOM ( 10 * 20 ) »GC0MR(5>999) rGCOMl (5»999) »pFXl(5»50) » 

3 PFY1 (5»50) » KPTS(l)»XCOMN(10r50)»YCOMN(10r50)»KNRCl)»KDR(l) 
DOUBLE PRECISION G 

IMPLICIT REaL*8 ( A— F • P-W ) 

REAL*9 X»Y» XV » CNR »CNI» XCOMN »YCOMN.CDR»CDI 
DIMENSION PPT (4) » LSN(4) 

COMMON TYPE (50) 

INTEGER TYPE 

COMPLEX PtPPT 

DO 140 J=l* NFREQ 

IF (TYPE (J) .EQ.»G')P=-PPT(1) 

IF (TYPE (J) •EQ. , P , )P=-PPT(2) 

IF(TYPE(J) .EQ. ’S’ >P=-PPT(3) 

IF (TYPE (J) .EQ. *A» )P=-PPT<4) 

IF(TYPE(J).EQ. »G’)SGN= LSN(l) 

IF(TYPE(J).EQ.»P')SGN= LSN(2) 

IF ( TYPE (J) »EQ. *S’ )SGN= LSN(3) 

IF(TYPE(J).EQ.*A»)SGN= LSN (4) 

KWHICH- KPTS(J) 

XV= OMEGA (KwHICH) * 6.2831853 
DO 130 L=l* KCHNL 
NCOMO= NUMC(L) 

KNOT=0 
LNOT=0 
IOP= KONT(L) 
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DO 130 N=1»NC0MD 
IF(N.GT.l)l0P=2 
Nl= NRATOR(L»N> + 1 
Ml= NDENOM(L#N) + 1 
DO 5 LP=1»N1 

5 X(LP)= XCOMN ( L r LP+KNOT ) 

DO 6 LP=lrMl 

6 Y(LP>= YCOMN(L»LP+LNOT> 

K2=N1-1 

K3=M1-1 

CALL P0LFV<X*K2>XV»CNRiCNI) 

CALL P0LFV(Y»K3»XV»CDR»CDI) 

RD= CNR**2 ♦ CNI**2 

RR= (CDR*CNR+CDI*CNI)/RD 

RI=(-CDR*CNI+CDI*CNR)/RD 

GR(J)= GC0MR(L#KWHICH)*RR - GCOMI (Lr KWHICH) *RI 
GI(J)= 6C0MR(LfKWHICH)*RI + GCOMI (LrKWHlCH) *RR 
A= REAL ( P ) +GOR ( KWH I CH ) -GCOMR ( L r KWHICH ) 

B- AIMAG(P)+GOI(KWHXCH)-GCOMI (LrKWHlCH) 

FREQ=1.0 

KSKIP=1 

DO 40 I=lrNl 

KUL1=(-1)**< (I + D/2) 

KUl2=(-1)**((I+2>/2) 

fuli=kuli 

FUL2=KUL2 

IF (KSKIP-1 )20»20r30 
20 C(1)=-GR(J)*FREG*FUL2 
0 ( I )=-GI ( J) *FREQ*FULl 
KSKIP=2 
GO TO 40 

30 C(I)=-GI(J)*FRE0*FUL2 
D ( I ) =-GR ( J ) *FREG*FUL 1 
KSKIP=1 

40 FREQ= FREQ*0MEGA<KwHICH)* 6. 2831853 
FREQ= 1.0 
DO 50 I=lrMl 
KMUL=(-D**( (I + l)/2) 
emul=kmul 

E(I)= -FREQ * EMUL 

50 FREQS FREQ * OMEGA (KWHICH ) *6. 2831853 
FNA1=0.0 
FNA2=0.0 
DO 60 IslrNl 
FNA1=FNA1+C(I)*XU) 

60 FNA2=FNA2+D(I)*X(I) 

FD2=0.0 

KI= 2 * UK3+D/2) 

DO 70 I=2rKlr2 
70 FD2=FD2+E(I)*Y(I) 

F01=0.0 

KE= 2 * ( (K3+2)/2) - 1 
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00 80 I=1*KE»2 
80 F01=FDl+E(I)*Y(I) 

FN1= FNA1 + FD1 * A - FD2 * B 

'FN2= FNA2 + FD2 ♦ A + FOl * B 

Fb=FDl**2+FD2**2 
FN=FN1**2 +Fn2**2 

FY£= (FD *<A * FNi + B * FN2> - FN * FDD/ FD**2 

FYO= (Fo*(-B * FNI + A * FN2> - FN * FD2)/ FD**2 

FX1=FN1/FD 

FX2=FN2/F0 

PFX1(L»KN0T+1)= 0.0 

DO 90 I=1»KE#2 

PFY1(L»I+LN0T)= FYE * Ed) * SON 
90 CONTINUE 

DO 100 I=2»KI»2 

PFY1 (L » I+LNOT ) = FYO * Ed) * SON 
100 CONTINUE 

IF(I0P.EQ.2)PFY1(L»LN0T+1)= 0.0 
DO 110 I=2*N1 

PFX1CL»I+KN0T)=(FX1*C(I) ♦ FX2*Dd>) * SGN 
110 CONTINUE 

KNOT= KNOT + N1 
LNOT= LNOT + Ml 
130 CONTINUE 
KLaO=0 

DO 135 IX=1»KCHNL 
KN0T= KNR(IX) 

DO 135 LX=l.KNOT 
KLAD=KLAD+1 

135 G(J.KLAD)= PFXKIX.LX) 

DO 139 IX=1»KCHNL 
LNOT= KDRdX) 

DO 139 LX=l»LNOT 
klad=klad+i 

139 G(J,KLAD)=PFYlUXfLX) 

140 CONTINUE 
NPaRC=KLAD 

DO 150 U=1*NFREQ 
SUM=0.0 

00 145 1=1»NPARC 
145 SUM=SUM+G(U»I)*G(J*I) 

SUM= DSQRT(SUM) 

DO 149 I=1»NPARC 

149 G(U»I)= G(J»I)/SUM 

150 CONTINUE 
RETURN 
END 


SUBROUTINE PHASEM ( GR » G I » KPO INT r NM r FQM IN , FQMaX » KPTS » STBM , OMEGA ) 
DIMENSION GR(1) »GI(1) #KPTS(1) »STBM(1) rOMEGA(l) 
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C SUBPROGRAM FOR CALCULATING PHASE MARGINS 
C 

C DEFINITIONS OF I/O VARIABLES 
C 

C GR -ARRAY OF OPEN LOOP TRANSFER FUNCTION REAL PARTS 

C GI -ARRAY OF OPEN LOOP TRANSFER FUNCTION IMAGINARY PARTS 

C KPCINT-NO. OF POINTS 
C OMEGA -ARRAY OF FREQS. 

C NM -COUNTER 

C KPTS -FREQUENCY NOS. WHERE MARGINS OCCUR 

C STBM -STABILITY MARGINS OF MARGINS 

C FQMIN -LOWER FREQ. FOR MARGIN DETECTION 
C FQMAX - UPPER FREQ. FOR MARGIN DETECTION 
P=I.O 

DO 50 I=1»KP0INT 
SO= GR(I)**2 + GI ( I ) **2 
S2=S0-1.0 
IF ( I .EG. 1 ) Sl=S2 
IF(OMEGA<I>.LT.FQMIN)GO TO 40 
IF (OMEGA ( I ) .GT .FQMAX ) RETURN 
IF l ABS (S2 ) .LT.1.0E-20) GO TO 30 
, SGN=S2/ABS<S2) 

IF (S1*SGN.GT .0.0) GO TO 40 
30 11=1-1 

IF(ABS(S2) .LT.ABS(Sl) )I1=I 
nm=nm+i 

KP1S(NM)=I1 

S3= (P+GR(H) >**2+GI (Il)**2 
stdm(nm)= sort (S3) 

40 S1=S2 
50 CONTINUE 
RETURN 
END 


SUBROUTINE GAINMG (GR»GI »KPOINT » NM * FQMIN » FQMaX » KPTS? STBM r OMEGA) 
DIMENSION GR(1) »GI ( 1 ) *KpTS ( 1 ) * STBMC 1 ) » OMEGA ( 1) 

C 

C SUBPROGRAM FOR CALCULATING GAIN MARGINS 
C 

C DEFINITIONS OF I/O variables 

c 

C GR -ARRAY OF OPEN LOOP TRANSFER FUNCTION REAL PARTS 

C GI -ARRAY OF OPEN LOOP TRANSFER FUNCTION IMAGINARY PARTS 

C KPOINT-NO. OF POINTS 
C OMEGA -ARRAY OF FREQS. 

C NM -COUNTER 

C KPTS -FREQUENCY NOS. WHERE MARGINS OCCUR 

C STBM -STABILITY MARGINS OF MARGINS 

C FQMIN -LOWER FREQ. FOR MARGIN DETECTION 
C FQMAX - UPPER FREQ. FOR MARGIN DETECTION 
P=1.0 
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DO 50 1 = 1 » KPOINT 
S2=GI(I) 

IF(I,EQ.l)Sl=S2 
IF(OMEGaU) .LT.FQMIN)GO TO 40 
IF (OMEGA < i ) .gt .fgmax)return 
IF(ABS<S2> .LT.1.0E-20)GO TO 30 
SGn=S2/ABS<S2> 

IF (S1*SGN • GT • 0 • 0 ) Go TO 40 
30 CR= GR ( I ) 

IF(CR.G£.0.0)G0 To 40 
11 = 1-1 

IF(ABS(S2).LT.ABS(S1) > 11=1 
nm=nm+i 

KPT S (NM) =11 

S3= (P+GR(U) )**2+GI(Il)**2 
STbM(NM)= SQRT(S3) 

40 S1=S2 
50 CONTINUE 
RETURN 
END 


SUBROUTINE SRMINS(GR»GI . KPOINT. NM.P .N.FQMIN.KPTS.STBM. OMEGA) 
DIMENSION GR(1) .GI (1) .KPTSU) .STBMU) » OMEGA (1.) 

SUBPROGRAM FOR DETERMINING THE MINMUNS OF THE OPEN LOOP FREQUENCY 
RESPONSE WITH RESPECT To THE -1 POINT, WHEN GIVEN POINTS ON THE 
OPEN LOOP REQUENCY RESPONSE 

DESCRIPTION OF I/O VARIABLES ' ,V 

GR - VECTOR OF REAL PARTS OF OPEN LOOP FREQUENCY RESPONSE 

GI - VECTOR OF IMAGINARY PARTS OF OPEN LOOP FREQ. RESPONSE 

KPOINT - NUMBER POINTS oF THE OPEN LOOP FREQ. RESPONSE GIVEN 
OMEGA - CORRESPONDING FREQUENCIES OF CHOSEN POINTS 
. KPT.S -FREQUENCY NOS. WHERE MARGINS OCCUR 
FQ'MIN -MINIMUM FRQ. CONSIDERED 
-P -POINT W.R.T. A MAX. OR MIN. IS DESIRED 
N -DETERMINES WHETHER A MAX. OR MIN. IS DETERMINED 

AS(n1 = 0.0 
S1=0.0 

DO 50 1=1. KPOINT 

IF (OMEGA ( I ) .LE.FQMIN) GO TO 50 

S2= (P + GR ( I ) ) **2 + GI(I)**2. 

ASN2=S2-S1 

5 CONTINUE . 

IF(ASN2*N)10. 50.10 
10 IF(ASNl*ASN2)20.40,40 
20 IF ( ASNl*N )30»40»40 
30 NM=j\|M+l 

u= i - i 

KPTS(NM)=I1 
STbM(NM)= SqRT(SI) 
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noooooooo 


40 S1=S2 

ASMs ASN2 
50 CONTINUE 
RETURN 
ENQ 


SUBROUTINE DIRVEC(G»NM»KPARC*DV* WEIGHT) 

directional VECTOR program 

DEFINITIONS OF I/O VARIABLES 

G -MATRIX WHOSE ROWS CONTAIN THE GRADIENT VECTORS OF THOSE 

stability margins only considered pertinent 

NM -NUMBER OF STABILITY MARGINS CONSIDERED PERTINENT 
WEIGHT-WEIGHTING FACTOR VECTOR 

DIMENSION G(20»99)> A(30>30)» WEIGHT ( 1 ) » 

/ AI(30.30)» X(30)» DV (30 ) » Y(30) 

IMPLICIT REAL*8 ( A-H t O-Z ) 

DO 200 K=1*NM 
Y(K)= WEIGHT(K) 

DO 200 J=K»NM 
SUMS 0.0 

DO 150 I— 1 » KPARC 
150 SUM= 'SUM '+ G ( J r I ) * G(KiI) 

A(j,K)= SUM 
A(K#J)= SUM 
200 CONTINUE 

IF(NM.GT.l>GO TO 230 
AI(1*1)= 1.0/A(lfl) 

X(l)= WEIGHT (1) * AI<1»1> 

GO TO 310 
230 CONTINUE 

CALL MATlNV(A»NMf AI»IER) 

IFllER.EQ.O)GO TO 300 
WRITE (6 » 250 ) 

250 format (*o**i5X»*the partials are not linearly independent, 
/he program is terminated.*) 
stop 

300 CALL MATMUL(NM»AI »NM»Y»1»X) 

310 CONTINUE 

DO 450 I=1*KPARC 

sum= o.o 

DO 400 J=1»NM 

400 SUMS SUM + G(J» I) * XU) 

450 DV ( I ) s SUM 
690 RETURN 

end 


SUBROUTINE MATINV(2,N»Y*IER) 


THUS T 
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00 50 1=1»KP0INT 
S2=GIU) 

IF ( I , EQ. 1 ) Sl=S2 
IF(QMEGA<I>.LT.FQMIN>GO TO 40 

IF(OMEGa(I> .gt.fqmax)return 
IF(ABS<S2) .LT.1.0E-20)GO TO 30 
SGN=S2/ABS(S2) 

IF(S1*SGN.GT .0.0) GO TO 40 
30 CR= GR<I) 

IFICR.GE.O.OIGO TO 40 
11 = 1-1 

IF(ABS(S2).LT.AaS(Sl) >11=1 

NM=NM+1 

KPTS(NM) = H 

S3= IP+GR(U))**2+GI(I1)**2 
STbM(NM)= SQRT (S3) 

40 S1=S2 
50 CONTINUE 
RETURN 
END 


SUBROUTINE SRM I NS ( GR # G I » KPO INT . NM » P * N * FQM IN , KPTS t STBM t OMEGA ) 
DIMENSION GR(1) »GIllJ *KPTS(1J »STBM(D »0MEGA(1> 

C 

C SUBPROGRAM FOR DETERMINING THE MINMUNS OF THE OPEN LOOP FREQUENCY 
C RESPONSE WITH RESPECT To THE -1 POINT WHEN GIVEN POINTS ON THE 
C OPEN LOOP REQUENCY RESPONSE 
C 

C DESCRIPTION OF I/O VARIABLES 

C GR VECTOR OF REAL PARTS OF OPEN LOOP FREQUENCY RESPONSE 

C GI VECTOR OF IMAGINARY PARTS OF OPEN LOOP FREQ* RESPONSE 

C KPOINT - NUMBER POINTS OF THE OPEN LOOP FREQ* RESPONSE GIVEN 

C OMEGA - CORRESPONDING FREQUENCIES OF CHOSEN POINTS 

C KPTS -FREQUENCY NOS. WHERE MARGINS OCCUR 
C FQMIN -MINIMUM FRQ. CONSIDERED 
C -p -POINT W.R.T. A MAX* OR MIN. IS DESIREO 
C N -DETERMINES WHETHER A MAX. OR MIN* IS DETERMINED 

ASNl=Q.O 
S1=0.Q 

DO 50 I=1*KP0INT 
IF(OMEGAU) .LE.FQMINJGO TO 50 
S2= (P + GRU)>**2 + 5I(I)**2 
ASN2=S2-S1 
5 CONTINUE 

IF l ASN2*N) 10» 50 * 10 
10 IF(ASN1*ASN2)20*40.40 
20 IF! ASNl*N) 30 » 40 >40 
30 NM=NM+1 
■11= I - 1 
KPTS(NM) = U 
STbM(NM)= SORT (SI) 
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noon on 


MULTIPLIES (A) * (B) 

A IS AN NR X N 
B IS AN N X NC 
X IS AN NR X NC 

00 4 I=1»NR 

4 X(I) = 0.0 
00 5 I=1»NR 
DO 5 K=1»NC 
DO 5 J=1»N 

5 X(I)= X(I) + A ( I . J) * B ( J) 
RETURN 

end 


SUBROUTINE POLFV ( Fw » K * X » FREAL .FIMAG) 

C PROGRAM FOR EVALUATING A POLYNOMIAL AT AN IMAGINARY FREQUENCY 
C 

C DEFINITIONS OF I/O VARIABLES 
C 

C FW -VECTOR POLYNOMIAL COEFFICIENTS 

C K -ORDER OF POLYNOMIAL 1 •' 

C X -FREQUENCY TO BE EVALUATED AT 
C FREAL -REAL PART OF Fw(JX) 

C FIMAG -IMAGINARY OF FW(UX) 
dimension Fw(1) 

KEVEN=<K+2>/2 

K0D0=(K+l)/2 

Y=I.O 

FREAL=0.0 

DO 10 I=1.KEVEN 

L=2*I-1 

FREAL= FREAL + FW(L)*Y 
10 Y=-Y*X*X 
FlMAG=0.0 

IF (KODD.EQ.O)GO TO 30 
Y=X 

DO 20 I=l»KODD 
L=2*I 

FIMAG= FIMAG + FW (L ) *Y 
20 Y=-Y*X*X 
30 RETURN 
END 


SUBROUTINE OTPT 1 ( STBM . OmEGA . KPTS » NM . XCOMN . YcOMN . KM IN » RQ . LOX . KCHNL . 
1 NUMC .NRATOR.NDENOM.PRX.PRY.il. 12) 

DIMENSION STBM ( 1 ) .OMEGA ( 1 ) .KPTS ( 1) » XCOMN ( 10. 50) » YCOMN { 10 .50) .RQ(1 ) 
1 »PRX(1) .PRY(l) .NUMC(l) »NRATOR(10.20) .NDENOM(10»20) 

DIMENSION TYPE (50) 

COMMON TYPE 
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INTEGER TYPE 
C 

C SUBPROGRAM FOR THE OUTPUT OF INFORMATION CALCULATED 
C 

IOP=l 

WRITE (6 >10) LOX 

10 FORMAT (* O' »25X» ’ITERATION NO. »>I4) 

DO 110 I=1*NM 
KWH= KPTS(I) 

FREO= OMEGA (KWH) 

IF(I.EG.KMIN+1)60 TO 50 
IF(I.EQ.l)GO TO 70 
GO TO 90 
50 WRlTE(6>60) 

60 FORMAT C O' >25X> 'ATTENUATED FREQUENCY INFORMATION 1 //) 

GO TO 90 
70 WRITE (6 >80) 

80 FORMATl *0' >25X> 'RELATIVE STABILITY INFORMATION'//) 

90 CONTINUE 

WRITE (6 > 100 ) I. STEM ( I ) > FREQ> RQ(I)> TyPEd) 

100 FORMATC ' »2X» 'MARGIN RADIUS NO, '>I2>»='>F10.5»5X>'FREQUENCY=*» 

1 F10. 5>1X, 'HZ' >5X> 'DESIRED MARGINS' >Fl0.5>5X> 'MARGIN TYPE=»*1X» 

2 Al) 

110 CONTINUE 

DO 104 I=1>KCHNL 
WRITE (6> 102) I 

102 FORMAT( »0» »25X> 'CHANNEL N0.'>12>' COMPENSATORS') 

L=NUMC(I) 

LX=1 

LY=1 

KX=0 

KY=0 

DO 104 IX=1>L 

KX=KX + NRATOR(I>IX) + 1 

KY=KY + NUENOM(I>IX) + 1 

WRITE(6> 106) (XCOMN(I>N) >N=LX>KX) 

WRI TE ( 6> 107 ) ( YCOMN ( I >N) > N=LY > KY ) 

106 FORMAT ( *0 ' >10X> 'NUMERATOR' >8F10.5) 

107 FORMAT < ♦ 0 » > 8X> 'DENOMINATOR ♦> 8F10 . 5 ) 

LX=KX+1 

LYsKY+1 
104 CONTINUE 

WRITE (6» 130 ) (PRX(I) rI=IOP»Il) 

WRITE (6> 120 ) (PRY ( I ) > I=IOP> 12 ) 

120 FORMAT ( ' 0 * » ' PART I Al_S WITH RESPECT TO Y».8El0.3) 

130 FORMAT* ' 0 '" PARTI ALS WITH RESPECT TO X»>8El0.3) 

WRITE(6>160) 

160 FORMATC O' ) 

RETURN 

END 
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SUbROUT INE MULOUT ( KCHNL , NUMC , NR ATOR » NDENOm , KNR ,KDR * XCOMN , YCOMN) 
DIMENSION Con (30 ) » COM(30), XCOF ( 30 ) » XCOmN (10,50) , YC0MN(10»50) , 
1 NUMC (30), NRATOR(10»20) » NDENOM ( 10 , 20 ) , KNR (20), KDR(20) 

DO 80 1=1, KCHNL 

Con ( l ) = l.o 
n=o 

LAX=1 

NAT- NUMC ( 1 ) 

WR1TE(6,40) I 

40 FORMAT ( »0* ,25X, 'CHANNEL NO. ’» 12, 2X, ’COMPENSATOR* ) 

DO 65 U=l, N at 
M=nRATOR(I»U) 

Ml= M + 1 
LAY= LAX + M 
KL=0 

DO 62 K=LAX,LAY 
KL=KL+1 

62 COM(KL)= XCOMN ( I ,K) 

LAX= LAY + 1 

CALL POLMULCCON, COM, N,M, XCOF) 

n=n+m 

N1=N+1 

DO 64 K=1,N1 

64 CON (K ) = XCOF (K) 

65 continue 

WR 1TE ( 6 , 67) 

67 FORMAT( 'O' ,25X, 'NUMERATOR COEFFICIENTS’) 

WRITE (6, 69) (CON(J) ,J=1,NI) 

69 FORMAT ( 'O’ ,2X,7E15.5> 

CON(1)=1.0 

N=0 

LAX=1 

DO 75 J=1 , NAT 
M= NDENOM ( I, J) 

Ml= M+l 
LAY= LAX+M 
KL=0 

DO 72 K=LAX,LAY 
KL=KL+1 

72 COM ( KL) = YCoMN ( I , K ) 

LAX= LAY + 1 

CALL POLMUHCON, COM, N»M, XCOF) 

N=N+M 

N1=N+1 

DO 74 K=l,Nl 

74 CON (K) = XCOF(K) 

75 CONTINUE 
WRITE(6,77) 

77 FORMAT('0»,25X, ’DENOMINATOR COEFIC IENTS* ) 

WRITE (6, 69) (CON(J) ,J=l,Ni) 

80 CONTINUE 
RETURN 
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ooooooooooooo 


'■ • END 


SUBROUTINE POLMUL ( CON » COM » N » M» XCOF ) 

DIMENSION CONU>» COM ( 1 ) » XC0F<1>» C0NA(5Q)# COMRA(50) 

FOR DOUBLE PRECISION REMOVE C FROM FIRST COLUMN OF NEXT CARD. 
DOUBLE' PRECISION CON# COM# XCOF# CONA# COMRA 

the vector con is a vector of the coefficient of a POLYNOMIAL 

OF ORDER N. 

the VECTOR COM IS a VECTOR OF THE COEFFICIENTS OF A POLYNOMIAL OF 
ORUER M. 

the vector XCOF is a vector of the coefficients of the product of 
A polynomial of order n and a polynomial of order m. the 

POLYNOMIAL of WHICH THE COEFFICIENTS ARE the VECTOR XCOF HAS AN 
ORDER of m + n. 

DO 1 1=1 # M 

1 CONA(i)=0.0 
NX=N+1 

DO 2 1=1 #NX 

LX=M+I 

2 CONA (LX ) =CON ( I ) 

MX=M+1 

DO 3 1=1 »MX 

MY=M+2-I 

3 COMRA(I)=COM(MY) 

DO 4 1=1. N 

NX=M+1+I 

4 COMRA(NX)=0.0 
KY=M+N+1 
KX=KY 

DO 7 K=1»KY 
XCOF(K)=0.0 
DO 5 L=1»KX 

5 XCOF (K)= CONA(L) * COMRA «L) +XCOF (K > 

KX=KX-1 

DO 6 J= 1 » KX 

6 CONA(vl)=CONA(J+l) 

7 CONTINUE 
RETURN 
END 
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APPENDIX B 


SUBPROGRAM SUMMARIES OF COMPENSATOR IMPROVEMENT PROGRAM 

INTRODUCTION 

In many instances , modifications of large computer-aided-design 
programs are necessary. This is especially true in cases where the 
program is to be adapted for solving problems other than those for which 
it was designed. In many situations these changes are to be made by 
someone other than the programmer who coded the program originally. 
Making changes in a program without prior knowledge of the theory and/or 
programming techniques used by the programmer can be a very time con- 
suming and laborious task. However, if certain specific and concise 
information is given, then, the amount of time and effort for changing 
or reproducing a program is decreased significantly. It is the purpose 
of this report to produce certain pertinent information concerning the 
subprograms of the CIP (Compensator Improvement Program) . With this 
information a programmer should be able to modify the programs or to 
replace any program with its equivalent. 
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Subroutine PARCLT 


The purpose of subroutine PARCLT is to calculate the gradients with 
respect to a control system's compensator coefficients of the distances 
squared between GH(jw) frequency response points and other points in the 
complex GH(jw) plane. The Vother points" are specified by the type of 
GH(joj) frequency response point under consideration, i.e., phase margin, 
attenuation margin, gain margin, or stability margin. Also, the points 
are chosen as pushing or pulling points according to their type. 

The equations for performing these calculations are given in the 
following. First consideration is focused on the general feedback con- 
trol system configuration shown in Figure 1. For this system with s = ju> 
the open loop frequency response is C(jw)/R(jw) when the feedback loop 
is broken at a. The compensated frequency response of the k^ channel 
when all channels are open is 

0 k 

— (ju) = e k (w) + • (!) 

til 

The k channel's compensator is assumed to consist of n^ sub-compensators 
in a cascaded arrangement or 

nk 

G (s) = n G,.,(s) (2) 

K i=1 Kl 


tlx til 

where G ki ( s ) * s th© * sub-compensator of the k channel. Each sub- 
compensator is assumed to have the following general form: 


V s > 


X_S C1 + X_ , S 11 " 1 + ... + X 0 


n 


n-r 


( 3 ) 


m 


m-1 


V + Vi s + ••• + y ° 
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Using the above notation the distance squared between some point, -P, in 
the complex GH(j«) plane and some point on the GH(jw) frequency response 
is represented as 


d(w) 


P + 


k=l 


[e k (»> 


Jf k («)] 


(4) 


This expression is now rewritten as a function of the p 1 "* 1 sub-compensator 
of the q t ^ 1 channel or 


d(oj) = 


A + jB + (c + jd)G qp (jto) 


(5) 


where 


A + jB = P + { l [e, (u>) + j f, (co) j } - 


k=l 


( 6 ) 


and 


{e (u) + jf (u>)} 

H 4 


c + jd = [e q (uO + jf q (u>)]/[G qp (jfi))] 


(7) 


Substituting (3) into (7) and manipulating results in 


n 


( X C i X i + A .l Q E 2j y 2j " B .l 0 E 2j+l y 2j+l } + 

k 


i=0 


n 


d(u>) 


( Jq D i X i + A Jo E 2j+l y 2j+l + *l Q E 2j y 2j } 


k 

I 

3=0 


( I E 2j y 2j ) + ( ^ 0 E 2j+l y 2j+l ) 


( 8 ) 


where k = m/2 and p = m/2-1 if m is even or k = (m-l)/2 and p = (m-l)/2 
if m is odd and the C’s and D's, and E's are defined by the following sets 


113 


[C 0 . 

Cl. 

Cz , 

C 3 , 

O 4 , 

C5, . 

••] = 

[c. 

-dto, -cto 2 , do> 3 . 

J- 

3 

O 

[D 0 , 

Pi, 

D 2 , 

P3, 

P 4 , 

D 5 , . 

..] = 

[d. 

co), -dw 2 , -cu 3 . 

du> 4 , 

[E 0 , 

El, 

e 2 . 

e 3 , 

E 4 , 

E 5 , . 

••] = 

[1, 

0 ), -u 2 , -w 3 , 01 4 

, W 5 , 


-dm 5 , . . . ] 

CW 5 ,...]: 

. . . ] . (8a,b,c) 


With the following definitions : 

n k 


rai " 


FN2 = 


J 0 D i X i + A J Q E 2 j+l y 2 j+l + B J Q E 2 j y 2 j 


FD1 = jFoVw 

F ° 2 = X E 2 j+l y 2 j+l 

J-u 

FD = (FDl ) 2 + (FD2 ) 2 
FN = „(FN1 ) 2 + (FN2 ) 2 


( 9 ) 


( 10 ) 

( 11 ) 


( 12 ) 

( 13 ) 

( 14 ) 


the partials of d(io) with respect to the q channel's p sub-compensator 

' - t * , 

coefficients are obtained as 


9d(w) 2 • [FD • (A • FN1 + B ‘ FN2) - FN \ FDl] • E e 

9 y e , .... .. . (FD ) 2 

for e even or . 


(15a) 


9d(u>) 2 • [FD • (-B • FN1 + A • FN2) - FN • FD2] • E 0 

S - - - - -- — 

9y g * • ; (FD ) 2 


(15b) 


for e odd and e = 1 , 2 , .... , m and 


9d(w) 2 • [FNl • C e + FN2 • D p ] 

•9x FN 

e 


( 16 ) 
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for e even or odd and e = 1,2, ...» n. Equations (8a,b,c), (9), (10), 

(11), (12), (13), (14), (15a, b) and (16) are programmed in PARCLT . The 
program has the capabilities of calculating the total gradient vector of 
the distance squared between some frequency response point and a chosen 
point in the complex GH(jw) -plane. 

The input and output variables and definitions for PARCLT is given 
in the following: 

Input Variables 

KCHNL - This is the number of channels fedback. It corresponds to 

j in Figure 1. 

NUMC(I) - This is a one dimensional array which specifies the number 

of sub-compensators in the I-th channel. 

NRAT0R(I , J) - A two dimensional array that denotes the numerator order of 

* 

the j-th sub-compensator in the I-th channel. 

NDEN0M(I,J) - A two dimensional array that denotes the denominator order 

* 

of the J-th sub-compensator in the I-th channel. 

XC0MN(I,J) - This is a two dimensional array which contains the numerator 
factors of the I-th channel’s compensator. The factors' co- 
efficients are listed in ascending order according to the 
powers of s where the s° coefficient is first and are listed 
succeedingly. The order and location of each factor is 
determined by NRAT0R(I,J). The J of XC0MN(I,J) denotes the 
J-th coefficient of all the numerator coefficients. The 
factor in which this coefficient belongs is determined by 

*NRAT0R(I,J) and NDEN0M(I,J) for a specified I arid J, which are equiva- 
lent, respectively, to q and p of (3), give the proper n and m of (3). 
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YCOMN(I,J) 


KNR(I) 


KDR(I) 


NFREQ 


KPTS(I) 


OMEGA (I) 


GOR(I) 


NRATOR(I,J) only. Once a particular factor and the location 
in XCOMN(I,J) is determined then the x's of (3) can be 
retrived from XCOMN(I,J). • 

- This is a two dimensional array which contains the 1-th 
channel's denomination factors in succeeding order. The 
factors are arranged in a parallel order to their orders, 
given in NDENOM(I,J). Using the orders in NDENOM(I,J) the 
location and length of a certain factor can be determined. 
Thus the y's of (3) can be retrieved from YCOMN(I,J) 

- A one dimensional array which contains the total number of 
numerator compensator coefficients of the factors in the 
I-th channel. 

- A one dimensional array which contains the total number 
of compensator coefficients of the denominator factors in 
the I-th channel. 

- The number of critical frequencies or frequency points for 
which gradients are to be found. 

- The frequency numbers of the NFREQ critical frequencies . 

If, for example, a frequency response is represented by 
348 frequency points then the frequencies are sequenced 
from 1 through 348. Thus, KPTS(I) contains the sequence 
number of the I-th critical frequency. 

- This one dimensional array contains the I-th frequency (in 
Hz) for representing the system. 

- A one dimensional array that contains the real part of the 
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GOI(I) 


GCOMR(I.J) 


GCOMI(I.J) 


KONT(I) 


open loop frequency response corresponding to the I-th 
* 

frequency. 

- A one dimensional array containing the imaginary part of the 

open loop frequency response corresponding to the 1-th 
* 

frequency. 

- The real part of the compensated frequency response of the 
I-th channel for the J-th frequency. This is equivalent to 
e i (ojj) in (1). 

- The imaginary part of the compensated frequency response of 
the I-th channel for the J-th frequency. This is equivalent 
to f ± ( w j) in (1)* 

- This one dimensional array specifies whether the d.c. gain of 
the compensator in the I-th channel is to be held constant. 

The d.c. coefficients (s° terms) of . the numerators of all sub- 
compensators are constrained to remain constant by automati- 
cally setting their partials to zero. This insures compensa- 
tor uniqueness. In every channel the d.c. coefficient partial 
derivatives of every sub-compensator except the first are 
automatically set to zero. Thus, the d. c. gain of each 
channel's compensator is assumed to be controlled by the d.c. 
term of the denominator of the first sub-compensator. If 

the partials for the I-th channel are being calculated and 
KONT(I) = 2, the partial of the d.c. term of the denominator 


Both GOR and GOI are related to the quantities in (4) by the following 
expression : 


j 

GOR(I) + jGOI(I) = l [e^) + jf k (« ± )] . 


(The j used to denoted S-T should not.be confused with summation termina- 
tion index, j). 
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TYPE(I) 


PPT(I) 


LSN(I) 


of the first sub-compensator of this channel is set to zero. 

If KONT(I) is something other than 2, this partial is 
unaltered. 

-A one dimensional array labeling the NFREQ frequencies. The 
quantities stored in TYPE are alphanumeric and consist of 
any of the letters G, P, S, or A which, respectively, stand 
for gain, phase, stability or attenuation margins. These 
s.ymbols are used to set the perturbation points (previously 
referred to as the "other points") and to set the sigh on 
the gradient. The gradient sign determines whether the per- 
turbation point is to be a pushing or a pulling point. 

- This is a complex one dimensional array that carries the four 
perturbation points respectively of gain, phase, stability, 
and attenuation margins (corresponds to P in Equation 4) . 

- This one dimensional integer array carries signs that deter- 
mines whether the perturbation points in PPT are to be pushing 
or pulling points. If LSN(I) = + 1 the point is a pushing 
point where if LSN(I) = - 1 the point is a pulling point. 


Output Variables 

G(I,J) - A real two dimensional array that contains the I-th critical 

frequency's scaled gradient vector with respect to the compen- 
sator's coefficients. The arrangement of every row of G is 
the scaled partials of all numerator coefficients starting 
with channel no. 1, sub-compensator no. 1 and progressing 
from sub-compensator to sub-compensator and from channel to 
channel until all scaled partials are listed. In the same row 
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of G following the numerator terms the scaled denominator 
partial derivatives are arranged similarly. The scale factor 
for the gradients is a quantity which converts each row of G 
to a unit vector. In a single tow of G there is an element 
for every compensator coefficient, even those whose partials 
are always set to zero. 

This variable specifies the total number of columns in G. 


Subroutine PHASEM 


The purpose of PHASEM is to detect and calculate phase margins of an 
open lpop control system which is represented by a discrete frequency re- 

i * 

sponse. The method for achieving this is given in the following discussion. 

It is assumed that the open loop frequency response is given in terms 
of real and imaginary parts. In particular suppose the i th frequency is f . 
then the corresponding real and imaginary parts of the open loop frequency 
response are GR^ and GI^. Phase margins occur at the real zeros crossing 
of the following sequence : 

2 

GR i + jGI i . (17) 

Next the following sequence is formed - 



Pi = S i * S i-1 * 


If <_'0 then or is zero or has made a zero crossing. Regard- 

less of which of these have occurred the phase margin frequency number is 
chosen as i or i - 1, depending on the smaller magnitude of or S^_^. 

The corresponding phase margin is calculated as S3 = 1.0 +GR^ + jGI^| 

where k is either i or i-1 as mentioned above. 

The input and output variables for this sub-program are as follows: 


Input Variables 

OMEGA(I) - A real one dimensional array that contains the frequencies in 
ascending order for describing the system. 

GR(I) - A real one dimensional array containing the I-th .real .part .of 
open loop frequency response. 
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GI(I) 

KPOINT 

FQMIN 

FQMAX 

NM 


NM 

KPTS(I) 

STBM(I) 


- This is a real one . dimensional array that contains the I~th 
imaginary. part of the open loop frequency response. 

- The number of frequency points used to describe the open loop 
frequency response of . the system* 

- The lowest frequency for which phase margins are to be detected. 

- The largest frequency for which phase margins are to be deter- 
mined . 

- NM + 1 is the number that the first margin found is to be 
given. For example, suppose NM is initially 2 and this program 
locates 3 margins. Then these margins would be labeled as 
margins 3, 4, and 5 respectively. 

Output Variables 

- This is the number that the last phase margin found is given. 

- A one dimensional Integer array that contains the frequency 
members of the margins found. 

- A one dimensional real array that contains margin values 
corresponding to the frequency numbers of KPT S. These margins 
are measured in terms of distances from the - 1.0 for jO.O 
point in the complex GH(ju) plane. 
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Subroutine GAINMG 


The purpose of this sub-program is to locate and calculate the gain 
margins of a discrete open loop frequency, response. The procedure used 
for accomplishing this is as follows. Suppose that the i C ^ frequency is 
f^. Then the corresponding real and imaginary parts of the open loop 
frequency response can be represented as GR^ and GI^ respectively. From 
the sequences of these real and imaginary parts the following sequence 
can be formed: 

U ± = GI ± • GI x (19) 

Whenever becomes negative or is zero a gain margin is detected. The 
frequency number of the gain margin is taken as i or i-1 depending whether 
|GI ± 1 > |GI ± _ 1 1 or |G1 ± j <_ |GI i _ 1 |. Then the gain margin is calculated as 

S3 = 

, <} 

For definitions of the I/O variables see Subroutine PHASEM. 


1.0 . + GR fc + 


JGI, 


where k is i or i - 1. 


Subroutine SRMINS 

The purpose of this sub-program is to calculate maxima or minima 
of a discrete open loop frequency response with respect to some chosen point 
along the real axis. Letting GR^ and GI^-be the 1 real and i 1 " imaginary 
parts, respectively, of an open loop frequency response, the following 
sequence is formed: 

2 

P + GR 1 + jGI i (20) 

where -P is some point along the real axis . From this another sequence is 
generated as follows: 



V i " U i- U i-1 


til 

If <_ 0 and > 0, the (i-1) frequency point corresponds to 

a relative maximum with respect to P. On the other hand, if <_ 0 

and ^ < 0 the (i-l) t ^ frequency point is a relative minimum with respect 
to P. 


The definitions of I/O variables are as follows: 


Input Variables 

GR(I) - See PHASEM. 

GI (I) - See PHASEM. 

OMEGA (I) - See PHASEM. 

KPOINT - See PHASEM. 


NM - See PHASEM. 

P - The negative of the real axis point for which maxima or 

minima are to be found. 
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N 


FQMIN 


- This integer variable determines whether the program is to be 
used for determining maxima or minima. Maxima are found if 

N = + 1, and minima are found if N = - 1. 

- The minimum frequency for which maxima or minima are determined 

* .v.' : 

All frequency points below this frequency are skipped. 

Output Variables 


Same as Subroutine PHASEM. 


-« Subroutine DIRVEC 

The purpose of this sub-program is to calculate the directional vector 
of the constraint improvement algorithm. The directional vector, d, is 

■* * V * \ ' 

calculated as 

d = [VG]a (22) 

where VG is a n x m matrix whose columns consist of the gradients of the 
active constraints. The quantity, a, is a m-componeht column vector 
which is determined from 

a = [VG T VG]" 1 c (23) 

where c is a m-component column vector whose elements are all positive. 
Definitions of I/O variables are given in the following lists: 

Input Variables 

G(I,J) - A two dimensional array whose rows are comprised of the 

gradients of the active constraints. 

NM - The number of rows in G. 

KPARC - The number of columns in G. 

WEIGHT (I) - A real one dimensional array that contains the column 

matrix c of (23). 

Output Variables 

DV(I) - A real one dimensional array which corresponds to d in (22). 
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Subroutine MATINV 

The purpose of this sub-program is to determine the inverse of a 

matrix. The method used is Gauss-Jordon reduction. It is assumed that 

— * 

no diagonal elements of the original matrix are zero. If in applying 
the Gauss-Jordon reduction procedure the magnitude of the i^ element of 
the i t ^ 1 pivot row is less than 1.0 x 10 -2 * it is assumed that the 
matrix does not possess an inverse. 

The I/O variables 'definitions are as follows: 

Input Variables 

A real two dimensional array whose inverse is desired. 

Number of rows and columns, in X. 

Output Variables 

A real two dimensional array that is the inverse of the X array. 

■ i • 

The error code of the program 
IER = 0 No error 

IER = 1 Matrix does not possess an inverse. 

1 ■ - : ! , ■ I. 


X(I,J) - 
N 


Y(I,J) - 
IER 
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Subroutine MATMUL 


The purpose of this sub-program is to multiply an n x,m matrix by 
a m x i column matrix. The equation for accomplishing this is 


m 

X i F % A ik b k 
1 k=l * 


(24) 


where A is the n x m matrix, b is the m-component vector (m x 1 column 
matrix), and x is the n-component vector resultant. 

The I/O variables for this sub-program are: 


Input Variables 

A(I , J) - A real two dimensional array (The matrix A in (24)). 

NR - An integer variable denoting the number of rows in A. 

N - An integer variable denoting the number of columns in A. 

B (I) - A real one dimensional array which contains the elements of b 

in (24). 

NC - Always chosen as the integer 1. 

X(I) - A real one dimensional array that contains the elements of 

x in (24). 
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Subroutine POLFV 


The purpose of this sub-program is to evaluate a polynomial at a 
point on the jw-axis of the complex s-plane. When given a polynomial 

n 

F (s) = J a. s (25) 

i=0 

where n is the order. The real and imaginary parts of F(s) when 
s = jo) are 

RE[F<Ju)] = l (-l)^ w 21 (26) 

i=0 

q 

IM[F(ju)] = l (-l) i a,,. 1 a) 2i+1 (27) 

i=0 

where p = n/2 and q = n/2 - 1 if n is even or p = (n-l)/2 and q = (n-l)/2 
if n is odd. 

Definitions of the I/O variables for POLFV are as follows: 

Input Variables 

FW(I) - A real one dimensional array that contains the coefficients of 

the polynomial that is to be evaluated. The arrangement of the 
coefficients is assumed to be in ascending order according to 
powers of s. 

K. - The order of the polynomial to be evaluated. 

X - The value along the imaginary axis for which evaluation is to be 

done [u in (26) and (27)]. 

Output Variables 

FREAL - RE[F (ju) ] 

FIMAG - IM[F(juO] 
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Subroutine 0TPT1 


The purpose of this sub-program is to output certain information at 
various stages of the main program. The information which is printed by 
this program is : 

1. The margin numbers 

2. The; frequency where each margin occurs 

3. The value of each margin 

4. The desired value of each margin 

5. The margin type, i.e., phase margin (P), gain margin (G) , 
stability margin (S), or attenuation margin (A) 

6. The directional vector at the last iteration 

7. The compensators at the last iteration. 
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Subroutine MULOUT 

The purpose of this sub-program is to convert the compensators which 
are in a cascaded factored arrangement into a single rational function 
form. It is assumed that the compensator for any channel is given by an 
equation such as 2. This compensator is converted to a single rational 
function form by multiplying all numerator factors together and multiply- 
ing all denominator factors together so that single polynomials are 
obtained for each. 

Definitions of the input and output variables for this program are 
the same as given in Subroutine PARCLT. 


Subroutine POLMUL 


The purpose of this program is to multiply two polynomials together. 
Assume there are two polynomials of the form: 


n , 
A(s) = l as 
i=0 

and 


B(s) = I b s 
i=0 


(28) 


(29) 


which are to be multiplied together. It is known that if (28) and (29) 
are multiplied together the resultant polynomial, P(s), will be of order 
m + n. Suppose that the coefficients of (28) are included as the last 
n + 1 elements of a vector which has m + n + 1 elements with the first m 
elements. as 0's. Denoting this vector as c, it becomes 


< C 1 


c , 
nr 


-m+1’ 


C m+n+l^ 


(30) 


where 


c^ = 0 i = l,2,...,m 


and 


c. = a. , i = in + 1; n + 2, .... i = n = 1 
i i-m-1 

Next , let the coefficients of (29) be cast in a vector d of the following 
form 


d 


Wj * ^2 * 


^mfl* ^nrt-2* 


^m+n+1^ 


(31) 
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where 


and 


d * “ b j . , i 1 & 2 % !••) m d* 1 

i m^i+1 * * 


d^ = 0 1 = m + 2, m + 3, ...,m+n+l 


(31) 


til 

From (30) and (31) the l coefficient of P(s) can be calculated as 


m+1 


p i ' . 1 “k+i 


(32) 


k=l 


i = 1, 2, . . . , m + n 


and 


m+n 

P(s) = l p,s . 
1=1 


(33) 


Definitions of I/O variables for subroutine POLMUL are: 

Input Variables 

CON(I) - A one dimensional array containing the coefficients of A(s) 

COM(I) - A one dimensional array containing the coefficients of B(s) 

N - The order of the polynomial, A(s). 

M - The order of the polynomial, B(s). 


Output Variables 

XCOF(I) - A one dimensional array containing the coefficients of the 
polynomial, P(s). 
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